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Chapter  I 

INTRODUCTION 

One  of  the  most  significant  contributions  to 
viscous  flow  theory  came  in  1904  when  Prandtl  [  1  ]  ^ 
published  his  boundary- layer  theory.  Since  that  time,  the 
theoretical  analysis  of  two-dimensional  boundary  layers 
has  been  advanced  to  the  point  where  accurate  solutions 
can  be  obtained  rather  routinely,  even  for  turbulent  flow. 
However,  most  flows  of  practical  interest  are  markedly 
turbulent  and  three-dimensional;  hence,  the  advancement  of 
three-dimensional  viscous  flow  calculation  methods  becomes 
important , 

An  effective  method  to  obtain  turbulent,  viscous 
flow  solutions  about  typically  encountered  aerodynamic 
configurations  is  to  numerically  solve  the  full  Navier- 
Stokes  equations  (using  some  means  to  model  the  turbulent 
Reynolds  stresses).  However,  this  approach  is  time- 
consuming  and  expensive  (sometimes  to  the  point  of  being 
prohibitive)  due  to  the  large  amount  of  computational 
resources  required  for  fully  three-dimensional  geometries. 


^Numbers  in  brackets  refer  to  similarly  numbered 
references  in  the  Bibliography. 
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As  an  alternative,  the  coupling  of  an  inviscid  flow  solver 
with  a  viscous  flow  (boundary-layer)  solver  has  proved  to 
be  a  useful  method  for  the  computation  of  viscous-inviscid 
interactive  flow,  particularly  two-dimensional  steady  flow 
(e.g.,  see  [ 2  —  4  J )  .  The  coupling  approach  has  also  been 
applied  successfully  to  three-dimensional  flow  [5].  The 
last  several  years  have  seen  significant  advances  in  the 
development  of  three-dimensional  Euler  equation  calcu¬ 
lation  methods  [6-9],  and  as  a  consequence  the  advancement 
of  three-dimensional  boundary- layer  solution  methods 
becomes  important  such  that  practical  three-dimensional 
problems  can  be  addressed  using  the  viscous-inviscid 
interaction  approach.  Therefore,  it  is  the  purpose  of  the 
present  study  to  develop  a  three-dimensional,  compressible, 
turbulent  boundary-layer  calculation  method  that  can  be 

used  for  transonic  flow  over  adiabatic  surfaces.  It  is 
emphasized  that  although  three-dimensional,  viscous- 
inviscid  interaction  is  the  ultimate  goal,  only  the  viscous 
(boundary- layer)  portion  of  a  coupling  approach  is  con¬ 
sidered  here. 

The  computation  of  three-dimensional  boundary 
layers  over  commonly  encountered  configurations  has 
received  considerable  attention  in  the  past  several  years 
(e.g.,  see  [10-20]).  Typically,  a  computational  method  is 
created  to  cater  to  a  particular  need  and/or  application. 
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The  following  discussion  addresses  the  reasoning  which  led 
to  the  present  approach. 

The  methods  described  in  [10-20]  can  be  classed  as 
either  integral  or  differential.  As  discussed  by  East 
[10],  and  Smith  [21],  among  others,  integral  methods  are 
computationally  faster  than  differential  methods  because 
the  former  have  one  less  space  dimension  to  contend  with 
and  also  because  more  empiricism  is  "built  into"  integral 
methods.  Although  generally  less  flexible  than  differ¬ 
ential  methods,  integral  methods  have  proved  to  be  as 
accurate  {and  in  some  cases,  more  accurate)  as  differen¬ 
tial  methods  for  two-dimensional  steady  flow  [22,23]. 
Therefore,  an  integral  approach  is  taken  here  in  the 
interest  of  speed,  and  also  because  the  three-dimensional 
method  is  an  extension  of  an  accurate  two-dimensional 
method  [24],  However,  the  question  of  accuracy  between 
integral  and  differential  methods  for  three-dimensional 
flow  has  been  investigated  to  a  much  lesser  degree  than 
two-dimensional  methods  [10,25,26]. 

When  using  a  coupled  approach,  it  is  desirable  that 
the  viscous  and  inviscid  surface  grids  (which  could  be 
nonorthogonal )  be  interchangeable  such  that  information 
generated  by  one  method  can  easily  be  conveyed  back  to  the 
other.  For  example,  surface  velocities  obtained  from  the 
inviscid  solution  must  be  used  as  input  to  the  boundary- 
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layer  equations.  The  steady  form  of  the  three-dimensional 
boundary-layer  equations  must  be  solved  on  a  grid  which  is 
dictated  by  domain-of-dependence  principles  [13,27] 
implying  that  interpolation  is  required  on  each  viscous/ 
inviscid  iteration.  However,  if  the  time-dependent 
boundary-layer  equations  are  used,  the  condition  of  inter¬ 
changeable  grids  can  be  achieved,  because  for  a  fixed  grid 
system  an  appropriate  time  step  can  be  chosen  to  maintain 
computational  stability.  Thus,  a  time-dependent  approach 
in  nonorthogonal  coordinates  is  adopted  in  order  to: 

(1)  provide  a  method  that  can  use  the  same  surface  grid  as 
an  inviscid  solver,  arid  (2)  account  for  time  accuracy,  if 
desired.  To  the  author's  knowledge,  all  previous  three- 
dimensional,  compressible,  turbulent  boundary-layer 
calculation  methods  have  been  for  steady  flow. 

The  system  of  equations  used  herein  is  the  three- 
dimensional,  time-dependent,  compressible  momentum  and 
mean-flow  kinetic  energy  integral  equations  in  non¬ 
orthogonal  curvilinear  coordinates.  To  take  advantage  of 
the  previous  work  in  two-dimensional  flow,  the  non¬ 
orthogonal  coordinates  are  related  to  streamline  coordi¬ 
nates  as  suggested  by  Smith  [14].  The  streamwise  velocity 
profile  used  is  that  of  Whitfield  et  al.  [28],  which  can 
be  expressed  analytically  over  the  entire  domain  0  £  < 


«e. 
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The  cross-flow  velocity  profile  used  is  the  triangular 
model  of  Johnston  [29]. 

In  this  study,  a  detailed  derivation  of  the 
boundary-layer  integral  equations  is  given.  The  necessary 
empirical  relationships  in  streamline  coordinates  are 
listed  in  addition  to  the  relations  between  streamwise  and 
nonorthogonal  quantities  for  the  Johnston  [29]  cross-flow 
velocity  profile.  The  numerical  scheme  used  and  stability 
and  convergence  for  various  spatial  difference  approxi¬ 
mations  are  discussed.  Finally,  computed  steady-state 
results  are  compared  with  measurements  and  with  compu¬ 
tations  of  previous  investigators. 
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Chapter  II 

DERIVATION  OF  EQUATIONS 

The  derivation  of  the  boundary-layer  integral 
equations  to  be  solved  is  given  in  this  chapter.  It  is 
important  to  include  the  derivation  of  these  equations 
because  they  apparently  do  not  exist  in  the  literature  for 
the  general  case  of  three-dimensional,  time-dependent, 
compressible  flow  in  nonorthogonal  curvilinear  coordi¬ 
nates.  An  abbreviated  version  is  given  in  [30],  but 
attention  here  is  focused  on  the  details  of  the  development. 


2.1.  Differential  to  Integral  Form 


The  differential  form  of  the  continuity  and 
momentum  boundary-layer  equations  for  steady  flow  is  given 
by  Cebeci  et  al .  [17],  These  equations  with  the  time- 
dependent  terms  are  given  below  in  a  compact  form: 


Continuity: 

h1h2  sinX 


M  +  _3_ 
at 


fpulh2  sinX)  + 
(pu3h1h2  sinX)  =  0 


(pu2h1  sinX) 

(2.1) 
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That  is,  when  i  =  1,  Eq.  (2.2)  is  the  x^-momentum  equation, 
and  when  i  =  2,  Eq.  (2.2)  is  the  X2-momentum  equation. 
Equations  (2.1)  and  (2.2)  are  written  for  a  general  non- 
orthogonal  curvilinear  coordinate  system  (fixed  in  time) 
like  that  depicted  in  Figure  1  (see  Nomenclature  for 
definition  of  terms).  It  should  be  noted  that  h3  =  1  in 


7 


7 


Figure  1.  General  Nonorthogonal  Curvilinear  Coordinate 
System  on  the  Body  Surface. 
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Eqs.  (2.1)  and  (2.2)  such  that  x^  is  the  actual  distance 
measured  normal  from  the  surface.  Also,  density,  pressure, 
and  velocity  appearing  in  Eqs.  (2.1)  and  (2.2)  are  time- 
averaged  turbulent  quantities,  and  t is  the  total  shear 
stress  (molecular  plus  turbulent)  in  the  x^-direction. 


2.1.1.  Momentum  integral  Equations 

Various  methods  can  be  used  to  arrive  at  the 
integral  form  of  the  equations.  The  approach  taken  here 
is  to  first  eliminate  the  pressure  gradient  terms  from 
Eq.  (2.2)  as  follows:  (1)  multiply  Eq.  (2.2)  by  sinX, 

(2)  write  this  result  at  the  edge  of  the  layer,  and  then 

(3)  subtract  this  result  from  that  of  Step  (1),  which  yields 

3u- _  3ua _  3ui 

Phih2  sinX  at-  +  cuih2  sinX  a3q;  +  pu2hi  sinx 

3“;  _  _2 

+  p  u3  hlh2  sinX  3^  +  phlh2  Ki+1  ui+l 


-2 


+  ph1h2  sinX  Ki/i+1  uxu2  -  phLh2  cosx  K.u. 

3u.  3u.  3ua 

-  ph1h2  sinX  ~  Puih2  ®inX  glT  "  pu2hl  SinX  3x2 


3u. 


-  pu3h1h2  sinX  ^  +  Phih2  cosX  Vi 

-  phlh2  Vl  Ui+1  "  phlh2  sinX  Ki,i+1  V2 

3tx. 

+  hlh2  sinX  vtf  =  0 


(2.4) 
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where  overbars  denote  boundary  layer  edge  values  and  the 
assumption  of  [TXj_]e£jge  -  0  has  been  made. 

Note  that  the  negative  of  the  eighth  through  the 
eleventh  terms  of  Eq.  (2.4)  can  be  written  as 

3u;  au_.  atu 

► 


d  i  dui  dui 

phih2  SinX  at-  +  pulh2  sinx  3X^  +  pu2hl  sinX 


3u  ■ 


+  pu3hih2  sinx 


*  ft  <Puihlh2  sinX)  +  3^  <Paluih2  sinX) 

+  W~2  lou2uihl  silU>  +  3^  (pUjUjhjhj  sinx) 


-  u 


3  3 

jjt  (Ph2h2  sin*>  +  (Pulh2  sin** 


+  3%  (pu2hl  slnXI  +  3^  <P“3hah2  sinX> 


(2.5) 


and  from  Eq.  (2.1),  the  term  in  brackets  in  Eq.  (2.5)  is 
zero  (metric  coefficients  are  assumed  to  be  invariant  with 
time).  Using  Eq.  (2.5)  in  (2.4)  and  rearranging,  the 
result  is 
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[p  hi 


3Ui 

h-  sin\  ‘ 


_  3ui  _  3u. 

an  +  ph2  sinx  ui  +  phi  sinX  u2  w: 


&u, 


+  pu3h1h2  Sinx  s^J-[ft  (Puihlh2  sinA) 

o 


0 

3*1 

(pu£u 

lh2 

9 

3x3 

(pu3u 

.  h.  h 

l  1 

hlh2 

Ki+1 

(pu 

hlh2 

sinx 

K. 

if 

Stx 

hlh2 

sinx 

ax. 

^2 


Q> 


=  o 


(2.6) 


<3> 


(The  numbered  brackets  appearing  in  Eq.  (2.6)  will  be  used 
later  on).  By  assuming  3u^/3x3  =  0,  adding  and  subtracting 
the  term 


3_ 

at 


(pUihih2 


sinX}  +  3^  (pUlUih2 


sinX ) 


+ 

in  Eq. 
in  Eq. 


3  g 

3*2  <  Pu2uih1  sinX)  +  a~5T^  (pu3uih1h2  sinx) 

(2.6),  and  using  Eq.  (2.1)  in  a  similar  fashion  as 
(2.5),  some  algebraic  manipulation  yields 
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9u, 


—  [hxh2  sinx  p(u.  -  u± )  ]  +  sinx(p-p)  ~ 


+  3^|h2  sinx  P“l(ui  -  V1  +  >hl  sinX  C^fUj-Uj)) 

_  au. 

+  h2  sinX(pu1  -  pux)  +  hl  sinx(pu2  -  Pu2) 


-  hlh2  cosX  K^pu*  -  pu?>  +  hjhj  Ki+1(pu?+1-pu2+1) 


+  sin;i1  Ki  i  +  l<pulu2  ~  pu^u2 ) 


+  nr  Ihih2 sinX  ou3(ui  -  V1 


a  fxi 

+  h.h„  sinX  - -  =  0 

12  9x, 


(2.7) 


Rearranging  the  first  line  of  Eq.  (2.7)  and  using  the 
identities 


— 2  2  —  —  — 

pui-pui  =  u.(pui-pui)+  pui(ui-ui)  (2.8a) 


pu 


lU2_pUlU2  =  ui(puj"Puj)+  Puj(ui-ui)  ,  i^j  (2.8b) 


results  in 
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h1h2  sinx 


3  /T77 


ft  ‘‘■V’V  -  ui  ft  <p-p> 

_  3u.  _ 

+  h2  sinX(pu^-pu^)  -^=-  +  sinx  (p  u2~pu2 }  ~ 

1  i 

+  tpulh2  »i»»<VV1  +  ^  [pujhj  aiMlIj-Ujll 

+  3x  tPu3hih2  sin^tui-ui)J  -  cosA  [ui (pui~pui ) 

+  pui(i.-u1)]  +  h^j  K.+1[5i+1(pui+1-pui+1) 


suj 


+  (mUl(ui+rui+l)1 


hih2  sinX  Kiii+i'ui(pui+i-pui+i1+  p'W'Wi 


3rXi 

+  h!h2  sinX  85TT  =  0 


(2.9) 


The  x^-  and  x2-moinentum  integral  equations  are  obtained  by 
integrating  Eq.  ( 2 . 9  >  over  0  <  x3  <  °°.  Using  the 
assumption  that  surface  metrics  are  independent  of  x^  [31], 
defining  the  following  integral  thicknesses  as 


00 


pqa K  =  ;  (pu,  -  pu.  )dx^ 

1  o 


pq^ij  y  puj(u,  -  u.)dx3 


(2.10a) 


(2.10b) 
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pq3Cij  =qS  pujlu?  -  u2)d*3 


(2.10c) 


q6* .  =  /  (u.  -  u. )dx, 
ui  n  1  1  J 


( 2 . lQd) 


pe  =  S  ( p  -  p  )dx. 


( 2.10e) 


and  considering  only  an  impermeable  wall;  i.e., 
*pu3*wall  =  °r  Eq*  becomes 


3 _ /  -  ~  <  _  ..  9 _ /  _  r 


pq2L' 


ui  ft<pep)J 

- —  [h7(p52  h2  sinA  eu) 

pq  h^h2  sinX  k-  1 


3  /Tt2 


+  hl  sinX  6i2} 


6*  3u .  6l  3u  •  (u. 

+  -k=  — -  +  -==  — -  -  K,  cotx[  =r-  6-  +  0.  ■ 

l  lij 


t^q  3x^  h2q  3x, 


+  Ki+i  c,cH"ir  6i+i +  ei+i.i+i 


*  Ki,i+ii=r  6l+i + 


-  i  c*xi  - 0 


(2.11) 


where  q  is  the  resultant  boundary-layer  edge  velocity  and 

Cf  is  the  local  skin  friction  coefficient  defined  as 
rxi 
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cfx 


i 


W 


(2.12) 


The  x^-  and  x2-momentum  integral  equations  result  from 
Eq.  (2.11)  by  setting  i  =  1  or  i  =  2,  respectively.  Again 
note  that  a  subscript  3  resulting  from  i  +  1  when  i  =  2  is 
taken  as  subscript  1.  Also,  a  comma  between  subscripts 
such  as  0^  is  taken  as  012  for  i  =  1  and  021  for  i  =  2. 

2.1.2.  Mean-Flow  Kinetic  Energy  Integral  Equation 

To  obtain  the  mean-flow  kinetic  energy  integral 
equation,  the  approach  is  to  multiply  Eq.  (2.6)  by  u^, 
integrate  both  over  0  £  and  then  sum  the  two 

resulting  integral  equations.  One  can  begin  this  program 
by  multiplying  Eq.  (2.6)  by  u.  and  writing  the  result  as 

u.(Lj  +  L2  +  L*  +  Lj)  =  0  (2.13) 

where  subscripts  1  to  4  in  Eq.  (2.13)  denote  the  terms 
included  within  the  similarly  numbered  brackets  of 
Eq.  (2.6),  and  subscript  i  on  u^  and  superscript  i  on  L1 
are  one  or  two.  Note  the  u^L*  term  of  the  above  can  be 
written  as 


UiL2 


■ui[ft{pui  hlh2  sinX)  +  3^(pului  h2  sinX) 

3  9  1 

+  3xJ(pu2ui  hl  sinX) + 3x^(Pu3ui  hah2  sinX)] 
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or,  using  the  continuity  equation. 


uiL2 


2[at(pui  hlh2  sinX3  +  ax^^l1^  h2  sinX) 

+  3^-i:Pu2ui  hi  sinX>  +  a^_{Pu3ui  hih2  sin^^2-14) 


Using  Eq.  (2.14)  in  Eq.  (2.13)  results  in 


_  8ui _  3u .  _  3u. 

Puih1h2  sinX  ^  +  pulUih2  sinX  ^  +  P^u^  sinX  ^ 


3t 


3u, 


+  pu3uih1h2  sinX  3x  2[3t(puihlh2  sinXJ 


+  8^(puluih2  sinX)  +  3l^(pu2uihl  sinX) 

g  2  |  _ 2  2 

+  ax“(  pu3uihlh2  sinX)J  ■  hih2  cosX  Kiui(Pui  "  Pui> 


+  h.h_  K.  .  u.(pu?  .  -  pu?  .  ) 
1  2  l+l  l  K  l+l  H  l+l 


+  hxh2  SinX  Ki>1+1  Ujlpu^j  -  pu1u2) 

9tx. 


+  hlh2  sinX  ui  9x- 


-  0 


(2.15) 


Adding  and  subtracting  the  term 


3  <  *77 2 


ft(Puihlh2  sinX)  +  3^(Puluih2  sinU 

+  3%*  Pu2Sihl  sinA)  +  8^'Pu3'Iihlh2  sinX) 
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to  Eq.  (2.15)  results  in 

,  _  au. 

pu;.h1h2  sin\  ^(pu^hg  sin**  +  puiuih2  sinX 

1  g  —0  _  ^ui 

•  I  a5^‘puiuih2  sin*>  +  pu2uihi  sinx  asq 

pu3uihih2  sinx  a5q 


\  3^(pu2uihl  s£nU 


2  3x3tpu3uihlh2  sin*^  +  2  {atlhlh2  sin£  plui  "  ui*l 

+  g^-[h2  sinX  pul(ui  “  “j1!  +  aJ-!*1!  sinX  pu2*“i  "  ui  *  J 
1  2 


+  sin*  pu3(u 


!-«?»} 


_ 2  2 

-  h1h2  cosX  Kiui(pui  -  pu^ 

+  hlh2  Ki+1  Vp“i+1  '  pui+l> 


+  sinx  Kijl+1  Ujlou^j 

3TXi 

+  h1h2  sinx  U1  -JJ-  =  0 


pulu2> 


(2.16) 


One  can  use  the  continuity  equation  and  rewrite  the  first 
eight  terms  of  Eq.  (2.16)  to  give 
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_  3^  _  _  9u . 

hlh2  sin^(Pu£  “  +  h2  sinA(pu1ui  -  pu.^)-^ 

_  8»i 

+  hl  sinX(pu2Ui  -  pu2u.)  — 

+  hxh2  sinX(pu3ui  -  P^^)^- 

+  l{ft^hlh2  sinAP(^i  “  ui>]+  sinApu^u?  -  u?)] 

+  3^[hl  *  ui)I+  d^[hlh2  •inApUjt^-u*)]} 

-  hjh2  oosx  W-,ul  -  puf>  +  hlh2  Ki+1  u.(FS?+1-pu^+1t 


pu.-pu.  =  u.(p-p)-  p(u.-u.)  (2.18a) 
p(u?-u?)  =  pui{ui-u.)+  u.lpu.-pu.)-  u?(p-p)  (2.18b) 
pujui-pujui  =  ^(pUj-pUj )-  puj(ui-ui)  (2.18c) 
ui(pui’Pui*  =  pui(u?-u?)+  u?(pui-pui  )-  pu^(u.-ui  )  (  2 . 18d  ) 
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ui ui+l~Pui+l 1 


,jui(ui+rui+il+  “i+i^vV 


-  PU1+1(W 


{ 2 . 18e) 


u  .  (n  u .u .-pu .  u  .  ) 
x  w  l  j  M  x  3 


pu..(u?-u?)  +  u?  (p  u  j-pu  j  ) 


-  pu.u.tu.^)  ,  i  7*  j 


( 2 . 18f  ) 


integrating  over  0  <  oo  using  the  integral  thickness 

definitions  given  in  Eq.  (2.10),  and  again  taking 

=  0,  results  in 


30^3X3  =  (pu,) 


3  wall 


Ui 


6uA  3u 

-f  0  - 
q3  p  q 


3  ”p-  -2k)jr  +  ^==3  n((”32  eii+  »  W 


'Vi  .  VuA  3U1  f  Vi  _  a2sh\  3uj 

vh1g2  h1q2  /  3X1  \h2q2  h2q2  /  3x2 


2h1h2  sinX  pq3LOAl 


Sx7<h2  slnX  P«3eil' 


-2 

u . 

1 


+  dblhl  sinX  Pq3  e12>]-  COtX  Ki  [£ii  +  =5  5I 


+  cscX  K 


i+1  1  i+1 , i  -2  "1 


-2  -2 

ui+l  P*  ui+l  r* 
+  - —  0  •  -  - 0„ 


-2  Uu; 


-2 

+  Ki,i+liEi,i+l  +  =1  6i+l 

SI 


u. u. . . 
1  l+l 


7? 


-2 

u . 

~i2  s 


u, 


1  r°°  9tx. 

I  u .  - — 

- 3  J  1  3x0 

P<3  n  3 


dx3  =  0 


(2.19) 
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The  mean  flow  kinetic  energy  integral  equation  is  obtained 
by  summing  Eq.  (2.19)  for  i  =  1  and  i  =  2,  yielding  the 
following  clean  but  formidable  equation: 
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Therefore,  in  summary:  (1)  the  momentum  integral  equations 
are  given  by  Eq.  (2.11)  for  i  =  1  and  i  =  2,  (2)  the  mean- 
flow  kinetic  energy  integral  equation  is  given  by 
Eq.  (2.20),  and  (3)  the  integral  lengths  are  given  by 
Eq.  (2.10).  All  metric  coefficients  are  listed  in 
Appendix  A. 

It  is  worth  noting  that  the  validity  of  the  deri¬ 
vation  just  presented  was  checked  for  some  special  cases 
by  comparing  Eqs .  (2.11)  and  (2.20)  to  corresponding 
equations  in  the  literature.  For  example,  Eq.  (2.11)  for 
i  =  1  and  Eq.  (2.20)  reduce  to  those  obtained  by  McDonald 
et  al .  [32]  for  the  case  of  two-dimensional,  time-dependent, 
compressible  flow.  Also,  Eqs.  (2.11)  and  (2.20)  reduce  to 
those  given  by  Nash  et  al.  [31,  Eqs.  (3.39),  (3.40),  and 
(3.44),  pp.  43  and  45]  for  the  case  of  three-dimensional, 
steady,  incompressible  flow  using  orthogonal  curvilinear 
coordinates  (i.e.,  A  =  tt/2).  However,  the  careful  reader 
will  note  that  a  sign  discrepancy  exists  in  one  term 
between  the  x^-momentum  integral  equation  as  given  by  Mager 
[33,  Eq.  (5.1),  p.  298]  and  Eq.  (2.11)  for  i  =  1  for  the 
case  of  three-dimensional,  time-dependent,  compressible 
flow  in  orthogonal  curvilinear  coordinates;  whereas, 

Eq.  (2.11)  for  i  =  2  and  Eq .  (5.2)  in  [33,  p.  298]  are 
identical  for  these  conditions.  After  careful  scrutiny, 
the  present  author  is  convinced  that  Eqs.  (2.11)  and  (2.20) 
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are  correct  and  that  the  sign  of  the  "curl^Q"  term  in 
Mager's  Eq.  (5.1)  [33,  p.  298]  should  be  negative.  The 
interested  reader  is  invited  to  assert  himself  concerning 
the  origin  of  this  discrepancy. 


2.1.3.  Restrictions  Pertaining  to  Steady  Edge  Conditions 
Equations  (2.11)  and  (2.20)  can  be  written 


y  po4i )-  ui  it 


i.  ,  i  *  1  or  2  (2.21) 


0— — 3  3 1 
2  p  q 


Pq2(eil  +  622  )+  pq(ul6I+U2S2)'  p6p(ul+u2  ) 


L  (2.22) 


where  and  L  are  defined  by  referring  to  Eqs.  (2,11)  and 

(2.20).  Expanding  the  derivatives  with  respect  to  time 
results  in  (for  i  =  1  and  i  =  2) 


36*  u.  36 

y  '  -  FE*  =  "  Tt. 

q  1 


(2.23a) 


u 


?t(0ll  +  022 


)  =  2qL  -Tl  (q«-1  -  ) 


1 
q 


u2  _ 

-T  ) 

q  *  *2 


u_  36 

=  3lT  =  **2 


2 


(2.23b) 


(2.23c) 
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where 


5-  .  0  un  a- 

1  3  /  7:771  p  1  i£ 

—  3t  (pq)  -  3t 

p  q  P  q 

3  .  %“2  3p 

p  q  p  q 


=  Z3{  1  611  +  922  ,Ie(p'>2)+5I  It  (pqV 


,  *  9  . - 

V  /  n  ^  t  t 


p  q 


+  6*  (P5u2)  -6p  ^  [p  (  u|  +  u|  )  ]} 


+  1 


fu 


Ae  !is  -,*0 

g  p  uJat  p  u^at 


(2.24a) 

(2.24b) 


(2.24c) 


Up  to  this  point,  Eqs.  (2.23a,b,c)  are  valid  for  time- 
varying  edge  conditions.  The  analysis  hereafter  is 
restricted  to  the  case  of  steady  edge  conditions;  that  is, 

M  -  33  -  lEi  -  !h  -  !He  _  „ 

3t  at  at  at  at 

and  therefore. 


This  is  a  physically  unrealistic  situation  in  that  the 
boundary-layer  edge  conditions  are  steady  but  the  boundary- 
layer  is  unsteady  developing  from  arbitrary  initial  con¬ 
ditions  that  do  not  correspond  to  reality.  However,  the 
solutions  were  found  to  be  insensitive  to  initial  con¬ 
ditions  as  discussed  in  Chapter  3,  Section  3.3.  With  these 
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restrictions,  Eqs.  (2.23a,b,c)  reduce  to 

3 5 1  u,  3 0 0  _ 

Tt  =“  3t”  *  q)ll  (2.25a) 

It*0!!  +  e22>  -  -  Ujli  -  V2  ( 2 . 25b  I 

_2__2__  =  q,2  (2.250 

Equations  (2.25a,b,c)  contain  the  13  integral 
lengths  6*,  0.  £■■,  5*  ,  and  6  where  i  =  1  or  2,  and 

-L  lj  1J  U-[  P 

j  =  1  or  2 .  In  the  present  analysis,  the  fundamental  step 
which  permits  the  construction  of  a  determinant  system  of 
equations  is  the  resolution  of  the  three-dimensional, 
turbulent  boundary-layer  velocity  profile  into  a  streamline 
coordinate  system  with  "streamwise"  and  "cross-flow"  com¬ 
ponents  as  illustrated  in  Figure  2.  The  streamline  coordi¬ 
nate  system  is  formed  by  the  projection  onto  the  body 
surface  of  the  external  streamlines  with  local  normals 
constructed  to  them;  the  direction  of  the  external  stream¬ 
line  is  called  the  "streamwise"  direction  and  the  "cross- 
flow"  direction  is  normal  to  it.  This  allows  each  velocity 
component  to  be  modeled  separately,  which  is  important 
because  it  has  been  observed  (see,  e.g.,  [34])  that  flow 
in  the  streamwise  direction  is  remarkably  similar  to  a 
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corresponding  two-dimensional  boundary  layer,  and  empirical 
relations  derived  for  two-dimensional  flow  (e.g.,  skin 
friction  and  shape  factor  correlations  and  velocity 
profile  families)  provide  good  approximations  to  the 
streamwise  components  of  velocity,  skin  friction,  and  so 
forth,  in  a  fully  three-dimensional  boundary  layer. 

Once  the  resolution  of  the  boundary-layer  velocity 
components  into  streamline  coordinates  has  been  accom¬ 
plished,  the  reduction  of  the  number  of  unknowns  appearing 
in  Eq.  {2.25)  is  hinged  upon  several  auxiliary  relations: 

(1)  empirical  relationships  in  streamline  coordinates, 

(2)  relations  between  streamwise  integral  lengths  using 
Johnston's  cross-flow  profile,  and  (3)  relationships 
between  streamwise  and  nonorthogonal  integral  lengths. 

The  following  sections  address  the  manner  in  which  these 
tasks  are  resolved  such  that  the  number  of  unknowns 
appearing  in  Eq.  (2.25)  is  reduced  to  three.  In  the 
following  discussion,  integral  lengths  written  with  upper¬ 
case  Greek  letters  represent  those  in  the  streamline  coor¬ 
dinates;  whereas  lower-case  letters  denote  integral  lengths 
resolved  in  the  nonorthogonal  system  (unfortunately,  this 
is  opposite  to  the  nomenclature  used  by  Smith  [14]). 
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2.2.  Empirical  Relationships  in  Streamline  Coordinates 

The  degree  of  success  of  an  integral  boundary-layer 
computational  procedure  is  ultimately  related  to  how  well 
the  auxiliary  relations  represent  reality.  The  present 
analysis  requires  models  for  (1)  the  streamwise  velocity 
profile,  (2)  the  cross-flow  velocity  profile,  (3)  skin 
friction  correlations,  and  {4)  shape  factor  correlations. 

2.2,1.  Streamwise  Velocity  Profile  and 

Shape  Factor  Correlations 

The  present  work  relies  heavily  upon  the  velocity 
profile  originally  postulated  by  Whitfield  [35]  for  steady, 
two-dimensional,  incompressible  or  compressible  adiabatic, 
turbulent  boundary  layers,  which  was  later  extended  to 
include  profiles  with  reversed  flow  [28,36].  It  should  be 
pointed  out  that  the  analytical  representations  of  the 
streamwise  and  cross-flow  velocity  profiles  do  not  appear 
explicitly  in  the  analysis;  rather,  shape  factor  corre¬ 
lations  and  relationships  between  the  streamwise  integral 
lengths  which  are  based  upon  the  velocity  profiles  are 
used. 

The  streamwise  velocity  profile  [28,36]  is  a 
function  of  "incompressible"  values  of  shape  factor,  H,  and 
momentum  thickness  Reynolds  number,  Ree  ,  that  is. 
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where 


(2.26) 


(2.27a) 


(As  seen  in  Eq.  (2,27a),  "incompressible"  here  means  simply 
that  integral  lengths  are  defined  to  be  independent  of 
density.  )  In  addition,  the  following  shape  factors  can  be 
defined  in  streamline  coordinates  as 


and 


0 


V  =  e 


li 


li 


0n 

V 


(2.27b) 

(2.27c) 

(2.27d) 


( 2.27e) 


where  the  above  streamwise  lengths  are  defined  in 
Appendix  B.  Using  Eq.  (2.26),  it  was  shown  in  [24]  that 
Hr  9i]/Qu'  and  H0*  can  be  correlated  with  H  and  edge  Mach 
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number,  Me,  with  only  a  weak  dependence  upon  Reg^.  The 
Hq*  correlation  actually  used  herein  is  based  upon  a 
streamwise  velocity  profile  valid  for  attached  and  sepa¬ 
rated  flow  [36].  All  shape  factor  correlations  used  are 
listed  in  Appendix  C,  including  the  Hq^  correlation 
originally  derived  by  Donegan  [37]  which  was  reported 
in  [30], 


2.2,2.  Cross-Flow  Velocity  Profile 

The  choice  of  cross-flow  velocity  profile  repre¬ 
sentation  can  have  a  significant  influence  on  the  results 
of  the  calculations,  as  shown  by  the  results  of  Smith’s 
integral  method  [14]  reported  by  East  [10],  Because  the 
main  objective  of  the  present  study  was  to  obtain  solutions 
of  the  three-dimensional,  turbulent,  integral  boundary- 
layer  equations,  it  was  decided  that  the  relatively 
uncomplicated  triangular  model  of  Johnston  [29]  would 
suffice,  with  the  understanding  that  an  improved  model 
could  be  used  later.  Smith  [14]  also  used  this  model  and 
obtained  reasonable  results  for  a  fairly  large  class  of 
three-dimensional,  turbulent  boundary- layer  test  cases 
[14,10,25,26].  However,  it  was  concluded  by  Johnston 
himself  [38],  that  "there  can  be  no  general,  universal 
cross-flow  model"  (for  example,  consider  a  "crossover"  or 
"s-shaped"  profile  where  the  cross-flow  velocity  is  both 
positive  and  negative  over  the  boundary-layer 
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thickness  [38]).  This  is  undoubtedly  one  of  the  weakest 
areas  concerning  the  use  of  integral  methods  for  the  calcu¬ 
lation  of  three-dimensional  turbulent  boundary  layers. 
Johnston's  cross-flow  model  [29]  is  given  as 


u 

n 


tan(ew) 


in  the  thin  layer  adjacent  to  the  wall,  and 


(2.31a) 


(2.31b) 


over  the  remaining  portion  of  the  boundary  layer,  where  A 
is  a  parameter  which  must  be  related  to  the  limiting  wall 
streamline  angle,  8W-  In  the  present  case,  the  relation¬ 
ship  originally  given  by  Johnston  [29]  and  later  modified 
by  Smith  [14]  as 


tan ( 8  )  =  A 
w 


0.1 


|c^  cos 


(bJ(  1+0.18  M=) 


TCI 


-  1' 


(2.32) 


is  used,  where  c^  in  the  above  is  the  local  skin  friction 
coefficient  resolved  along  the  external  streamline  flow 
direction  (in  the  present  case,  is  solved  for  itera¬ 
tively  by  knowing  A,  c^,  and  Mg) . 

2.2.3.  Skin  Friction 

The  skin  friction  correlation  used  is  that  given 
by  Whitfield  et  al.  [28],  as 
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= 


0. 3e 


-1.33H 


(  logio^en 


1. 74+0.31H 


(l.lxlO-4,  [ta„h(4  -  o^75  )-l] 


(2.33) 


where  c^,  H,  and  Reg^  denote  "incompressible”  values. 

The  first  term  on  the  right-hand  side  of  Eq.  (2.33)  was 
derived  by  White  [23]  (Eq.  6-179,  p.  518)  from  curve  fits 
of  the  Law-of-the-Wal 1/Law-of- the-Wake  using  Coles’ 
constants  in  the  Law-of-the-Wal 1  [23].  The  second  term  of 
Eq.  (2.33)  was  originally  reported  in  [39]  and  later  in 
[28]  and  [36]  and  was  appended  to  White’s  relation  to 
allow  c^  to  become  negative.  Although  separated  flows  are 
not  addressed  in  this  work,  Eq.  (2.33)  is  used  because  of 
its  behavior  at  high  shape  factors  (for  low  shape  factors, 
the  contribution  of  this  term  is  negligible) . 

The  relations  used  to  interrelate  compressible  and 
incompressible  variables  in  Eq.  (2.33)  are  Coles’  "Law  of 
Corresponding  Stations"  [40]  given  as 


Re.  cf  =  Ren  cf  (2.34a) 

0H  f  0  n  f 

and  the  correlation  of  Winter  and  Gaudet  [41],  which 
relates  c^  to  c^  by  the  relation 
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or,  in  effect 


where 


Re 


0 


11 


Re 


0 


11 


(2.34b) 


(2.34c) 


{ 2 . 34d  > 


Finally,  the  equations  used  to  resolve  cf  in 
streamwise  coordinates  into  nonorthogonal  components 
(which  Eqs.  (2.25a,b,c)  contain)  are  those  used  by 
Myring  [13]  and  Smith  [14], 


cf  =  cf 
X1 


sin(C)  -  cos(C)  tan  ( 3W) 
sin  ( X) 


(2.35a) 


cf  _  cf 
x2  1 


sin  (a)  +  cos  (a)  tan  (£w) 
sin  { X) 


(2.35b) 


where  Cf  is  the  value  of  skin  friction  resolved  in  the 
xi 

x^-direction,  a  is  the  angle  between  the  local  resultant 

edge  velocity  vector  and  the  x^-axis,  X  is  the  angle 
between  the  x^-  and  X2~axes,  and  £  =  X-a. 
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2.2.4.  Dissipation  Integrals 

As  shown  in  Appendix  D,  the  dissipation  integrals 
appearing  in  Eq.  (2.20)  can  be  written  in  terms  of  stream- 
wise  and  cross-flow  velocities  and  then  integrated  by 
parts  using  Johnston's  cross-flow  profile  [29],  yielding 


where  D^,  D^,  t^,  and  t2  are  defined  in  Appendix  D.  The 
"streamwise"  dissipation  is  evaluated  in  this  study 
using  the  correlation  developed  by  Donegan  [ 42]  and  later 
improved  by  Thomas  [ 43] ,  Actually,  the  product  c^  D^/2 
was  correlated  as 


c 


(2.37) 


This  correlation  was  derived  by  numerically  evaluating 
D®  using  a  constant  laminar  plus  turbulent  shear  stress  in 
the  region  very  near  the  wall,  the  Cebeci-Smith  eddy 
viscosity  model  [44]  in  the  inner  and  outer  regions,  and 
the  derivative  of  the  velocity  profile  used  in  [28].  In 
the  present  study,  the  contribution  of  has  been 
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neglected  in  comparison  to  with  the  understanding  that 
significant  error  could  be  introduced  in  flows  with  large 
crossflow  (see  Appendix  D). 

2.3.  Relationships  Between  Streamwise  Integral 
Lengths  Using  Johnston's 
Cross-Flow  Profile 

As  defined  by  Eq.  (B.8)  in  Appendix  B,  the  stream- 
wise  integral  length  a*  is  given  by 

CO 

pqA~  =  -  /  pu  dx,  (2.38) 

A  0  n  j 


Using  Eq.  (2.31b)  in  Eq .  (2.38)  (according  to  Smith  [14], 
only  the  outer  part  of  Johnston's  model  is  needed  to 
evaluate  streamwise  integral  lengths)  results  in 

00 

1>  q/\  =  -  A  /  (  pu  -  pu  )  dx- 

^  g  S  S  J 


-  -  A  /  [  (  p  u  -  pu  )  -  u  (  p  -  p  )  ]  dx~ 
o  *  ^  b  o 

=  -  A  [  p qA*  -  p us  ep  ] 

Thus, 

A*  =  -  A (AJ  -  ep  )  (2.39) 
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It  is  shown  in  Appendix  E  that  the  remaining  integral 
lengths  in  streamline  coordinates  (012'  021'  etc*^  can 
also  be  related  to  A,  A*,  A*,  01X,  En  ,  0U»  and  0^. 

2.4.  Relationships  Between  Streamwise  and 
Nonorthogonal  Integral  Lengths 

As  pointed  out  by  Smith  [14],  all  integral 
quantities  in  one  axis  system  are  uniquely  related  to 
those  in  the  other  and  these  relations  are  independent  of 
the  choice  of  cross-flow  profile.  For  example.  Smith  [14] 
shows  that  the  momentum  thickness  0^  in  the  nonorthogonal 
system  is  related  to  those  in  the  streamline  system  as 

0  =  — ^2“[011  sin2£-(Q12+02l)sin^  COS^+022  cos2Sl  (2.40) 

sin  X 

Stock  [15]  has  listed  all  of  these  relationships  using 
different  nomenclature  than  that  used  herein;  whereas 
Myring  [13]  and  Smith  [14]  give  the  relationships  between 
streamline  and  nonorthogonal  displacement  and  momentum 
thicknesses.  The  complete  list  using  the  present  nomen¬ 
clature  is  given  in  Appendix  F. 

2.5.  Reduction  of  the  Number  of  Unknowns--Summary 

The  number  of  parameters  in  the  system  can  now  be 
related  to  the  three  unknowns  9^,  H,  and  A.  This  process 
can  be  summarized  as  follows.  If  the  values  of 
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and  A  are  known  at  all  points  on  the  surface  and  at  a 
particular  time,  and  given  the  edge  conditions  q,  u^,  u^, 
Me»  and  geometric  parameters  a,  X,  and  £,  the  following 
sequence  of  calculations  produces  all  remaining  unknowns 
in  the  system  (Eq.  (2.25)): 


1.  Compute  shape  factor  correlations 

_  % 


H  =  H(H,Me) 


Hq  *  =  He*(H,Me) 


H0O  ■  vh'V 


0 


u 


0 


11  W11 


WT7  <H'Me> 


(see  Appendix  C) 


A1  =  H  011 


En  "  He*  0n 


2.  Compute  streamwise  integral  quantities 

021f  ®12r  etc*  (see  Appendix  E). 

3.  Compute  streamwise  skin  friction  (Eq.  (2.33)), 

Bw  (Eq.  (2.32)),  c^D^/2  (see  Appendix  C), 

Cf  and  Cf  (Eq.  (2.35)). 

X1  x  2 

4.  Compute  nonorthogonal  integral  quantities  6^, 
^2'  ®ii'  etc.  (see  Appendix  F). 
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2.6.  Formulation  for  Solution 

Using  the  relations  between  streamwise  and  non- 
orthogonal  integral  lengths  given  in  Appendix  F,  the 
aforementioned  empirical  correlations,  and  assuming  steady 
edge  conditions,  it  is  shown  in  Appendix  G  that 
Eqs.  (2.25a,b,c)  can  be  recast  into  matrix  form  as 

A  ^  (JJ)  =  b  (2.41) 

where 

A 

A  =  3x3  coefficient  matrix 
U  =  vector  of  unknowns 

*  (en,  h,  a)t 

b  =  RHS  vector  containing  spatial  derivatives 
and  edge  conditions 

The  elements  of  the  matrix  A  and  vector  b  are  also  given 
in  Appendix  G.  Equation  (2.41)  represents  a  system  of 
three  first-order,  nonlinear,  coupled  partial  differ¬ 
ential  equations  for  the  three  unknowns,  H,  and  A. 

The  numerical  approach  taken  here  is  to  first  reduce  the 
above  system  of  partial  differential  equations  (PDE)  to  a 
system  of  ordinary  differential  equations  (ODE)  and  then 
use  a  standard  integrator  for  ordinary  differential 
equations.  Particular  aspects  of  the  numerical  method 
used  herein  are  discussed  in  Chapter  III. 
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Chapter  III 

NUMERICAL  METHOD 

The  aforementioned  numerical  approach  is  commonly 
referred  to  as  the  Method  of  Lines,  which,  as  discussed  by 
Ames  [45,  pp.  302-304],  is  primarily  Russian  in  origin 
(see,  e.g.,  Liskovets  [46]).  This  numerical  method  was 
chosen  for  the  present  study  because  of  the  success 
demonstrated  by  Jameson  et  al.  [47]  who  used  this  approach 
to  solve  the  unsteady  Euler  equations  in  transonic  flow, 
and  in  addition,  showed  that  when  using  a  Runge-Kutta 
scheme,  convergence  to  a  steady-state  solution  could  be 
accelerated  by  using  a  local  time-step  as  dictated  by  the 
local  CFL  number  (Courant-Friedrichs-Lewy  stability 
criterion).  The  approach  taken  here  is  also  a  Runge-Kutta 
scheme  using  a  variable  time-step  to  accelerate  convergence 
(transient  results  were  not  considered  important  for  the 
cases  presented). 

3.1.  Implementing  Four-Stage  Runge-Kutta 
for  a  System  of  PDE's 

Consider  the  single  linear  model  equation 

|H+a|^=0,  a=  constant  >0  (3.1) 

d  l.  3  X 
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applied  in  the  x-t  plane.  By  discretizing  the  continuous 
x-t  space  as  depicted  in  Figure  3a,  Eq.  (3.1)  can  be  con¬ 
verted  to  an  ordinary  differential  equation  by  expressing 
the  spatial  derivative  with  an  appropriate  finite- 
difference  approximation,  or 

du . 

■qT—  =  -  aD  (u?)  =  f(t,u?)  (3.2) 

at  xi  l 

where  D  (u?)  is  a  finite-difference  operator  and  i  and  n 
denote  the  ith  value  of  u  at  a  time  level  n.  Writing 
Eq.  (3.2)  at  each  i  mesh  point  at  the  nth  time-level 
results  in  a  system  of  ordinary  differential  equations 
which  can  be  integrated  using  any  standard  ODE  integration 
scheme. 

As  mentioned  above,  a  Runge-Kutta  (R-K)  scheme  was 
chosen  for  the  present  study.  The  particular  scheme  is 
what  is  usually  referred  to  as  the  "classical"  R-K  scheme 
(Eq.  (20),  p.  120  of  Lambert  [48])  is  given  here  as 

u?+1  =  u”  +  ^  +  2^2  +  2k3  +  k4^  (3.3a) 

where 

k1  =  f(  t,u?  )  (3.3b) 

k2  =  f(  t  +  ^  ,  u?  +  (3.3c) 
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k3  =  f(  t  +  A|  r  u”  +  k2  )  (3.3d) 

k.  =  f(  t  +  At  ,  u?  +  At  k,  )  ( 3 . 3e ) 

4  l  j 


which  advances  the  solution  at  each  i  mesh  point  at  time- 
level  n  to  time-level  n+1  with  a  local  truncation  error  of 
0( At5) . 

For  a  model  equation  containing  two  space 
dimensions,  Eq.  (3.1)  becomes 


Su 
3 1 


+ 


a 


+  b 


3u_ 

3x2 


0 


(3.4) 


where  a  and  b  are  positive  constants.  Using  a  space-time 
discretization  as  depicted  in  Figure  3b,  the  corresponding 

ODE  is 


dt 


-  aD„  ( 


n 

Uij 


)  -  bD„  ( 


n 

u .  . 
13 


(3.5) 


n 

g(  truL. 


Applying  the  R-K  scheme  at  each  (i,j)  surface  grid  point 
results  in 


u 


n+1 
i  3 


(k^  +  2k2  +  2k^  +  k^  ) 


where  now. 


(3.6a) 
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kx  =  g(  t,u^  )  (3.6b) 
k2  =  g(  t  +  ^  ,  uj.  +  ^  k1  )  (3.6c) 
k3  =  g(  t  +  ,  u?_.  +  4|  k2  )  (3.6d) 
k4  =  g(  t  +  At,  uj  +  At  k3  )  (3.6e) 

The  present  set  of  equations  (Eq.  (2.41))  can  be  written 
in  vector  form  as 


(3.7) 


(3.8) 


M  and  N  are  3x3  matrices,  and  c  is  a  vector.  Thus, 
replacing  the  in  Eq.  (3.5)  with  results  in 


^~ii 


-  M. 


13  xn 


(  uV  )  -  N  -  -D  (  U‘.  .  )  +  c  (3.9) 

Xj  A «  J 


and  the  solution  for  IL  ^  can  be  advanced  from  n  to  n+1 
tising  Eq.  (3.6)  for  each  unknown  (i.e.,  0^  ,  H,  and  A). 
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Although  a  system  of  equations  is  involved,  the 
numerical  method  described  above  is  explicit  (i.e.,  all 
quantities  on  the  right-hand  side  of  Eq.  (3.6a)  are  known). 
For  stability,  explicit  methods  are  typically  restricted  • 
to  CFL  numbers  less  than  one.  However,  as  discussed  in 
the  next  section,  the  R-K  scheme  used  herein  permits  the 
solution  process  to  advance  using  CFL  numbers  greater  than 
one,  depending  upon  the  type  of  space  differences  used. 

3.2.  Stability  and  Convergence 
3.2.1.  General 

As  pointed  out  by  Mitchell  and  Griffiths  [49, 
pp.  181-182],  "a  stability  analysis  of  a  hyperbolic  system 
in  two-space  dimensions  is  extremely  difficult,  even  with 
A  and  B  constant"  (where  in  this  case,  A  and  B  correspond 
to  and  j ) .  For  example,  see  [50]  and  [51].  It 
follows  that  stability  requirements  used  in  the  present 
study  do  not  stem  from  an  analysis  on  the  complete  system 
of  equations  in  two-space  dimensions,  but  rather  from  a 
Fourier  analysis  based  upon  a  linearized  version  of 
Eq.  (3.9)  which  has  been  "split"  into  two,  one-dimensional 
(in  space)  problems.  That  is,  stability  requirements  are 
based  upon  the  "worst  case"  stemming  from  a  separate 
analysis  of  each  of  the  "split"  equations 
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and 


dU,  . 

=  -  M .  . D  {  U1?  .  ) 

at  13  x. 


dU.  . 

-^2-  =  -  N  .  .  D  (  Un  .  ) 
dt  i]  x„  -13 


(3.10a) 

(3.10b) 


Actually,  the  requirements  for  stability  come  from 
analyzing  the  single  linearized  two-dimensional  model 
equation,  Eq.  (3.5),  where  a  and  b  are  interpreted  as 
eigenvalues  of  the  and  matrices,  respectively 
(see  Lambert  [48],  p.  227).  Therefore,  the  equations  from 
which  the  stability  criterion  is  derived  reduce  to 


and 


du . 

_ 13  _ 

dt 


=  — p • • (M) D  (u?,) 


13 


X, 


13 


(3.11a) 


du .  . 

=  -Pij(N)Dx^(u”j)  (3.11b) 

where  p..(M)  and  p..(N)  are  the  spectral  radii  of  the  M.. 

x  J  1 3  1  j 

and  Ni_.  matrices,  respectively,  at  each  surface  grid  point. 
A  stable  time-step  is  chosen  as  follows:  at  the  beginning 
of  each  time-step  and  for  each  mesh  point: 

1.  Compute  the  elements  of  the  M  and  N  matrices. 

2.  Compute  all  eigenvalues  of  both  M  and  N 
matrices  and  determine  the  spectral  radius  of 
each;  that  is,  find  p^M)  and  pi  ^  ( N )  . 
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As  will  be  shown  presently,  a  local  CFL  stability  criterion 
for  the  four-stage  R-K  scheme  is  given  by 


p(M)  or  p(N) 


4*1  <  CFL 

Hij 


where  CFL  is  determined  from  a  stability  analysis  of  the 
R-K  scheme,  and  Ax  is  either  Ax^  or  Ax2«  Thus, 

3.  Compute  each  local  time-step  as 


Compute  At^j  =  min[ 


( AtXl) 


( AtX2) 


Although  not  mathematically  rigorous,  this  analysis 
for  determining  a  locally  stable  time-step  has  proved  to 
be  somewhat  conservative,  at  least  for  the  cases  con¬ 
sidered.  For  example,  numerical  experiments  revealed  that 
the  van  den  Berg  and  Elsenaar  test  case  (Section  5.4} 
could  be  computed  successfully  using  backward  space  differ¬ 
ences  and  a  local  time- step  which  would  require  a  CFL  of 
1.55  (it  will  be  seen  that  linear  stability  analysis  gives 
a  maximum  stable  CFL  of  approximately  1.3).  It  seems 
plausible,  however,  that  just  the  opposite  could  occur 
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when  considering  the  liberties  taken  during  the  course  of 
the  analysis;  i.e.,  the  "splitting"  of  the  complete  problem 
into  two,  one -dimensional  problems,  and,  of  course,  using 
a  CFL  number  derived  from  a  locally  linearized  equation. 
Therefore,  time-steps  as  determined  from  the  above  analysis 
should  be  used  with  care. 

3.2.2.  Stability  and  Convergence  Using 

Various  Space  Difference  Approximations 

The  stability  characteristics  of  the  R-K  scheme 
used  herein  depend  explicitly  upon  the  type  of  difference 
used  to  approximate  the  spatial  derivative  (Dx^  and  DX2 ) 
in  Eqs.  (3.11a,b).  Stability  of  the  model  problem 
Eq.  {3.2}  is  investigated  in  the  present  study  using  a 
Fourier  analysis  (45,  pp.  47-48]  with  three  different 
representations  of  D^fu?)  ((1)  first-order  backward, 

(2)  second-order  central,  and  (3)  second-order  backward).'*' 
Using  this  method,  the  propagation  of  a  single  row  of 
errors  represented  by  a  Fourier  series  along  some  initial 
time  line  is  examined  to  investigate  the  growth  of  the 
error  for  large  values  of  time. 


^In  the  present  work,  the  "order"  of  a  difference 
approximation  is  to  be  interpreted  as  the  order  of  the 
local  truncation  error.  For  example,  a  second-order 
difference  has  a  local  truncation  error  of  0(Ax2). 
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Using  a  first-order  backward  difference  for 
3u/3x  in  Bq.  (3.1)  results  in 


d  (U?)  ;  _i  (  un  -  u?  ,  t 

x  1'  ax  1  l-l 


(3.12) 


Thus,  Eq.  (3.2)  becomes 


H  V(u?)  =  f{t'u?> 


(3.13a> 


where 


V(u") 


n  n 

u  -  -  u .  , 
i  l-l 


(3.13b) 


Similar  relations  are  obtained  for  second-order  central 


D*<u?>  =  255  '  “HI  -  Ui-1  > 


(3.14) 


and  second-order  backward  differences: 


_  ,  n ,  •  1  ,  n  .  n  ,  n  . 

Vui>  ■  2Ei  '  3ui  -  4ui-l  +  ui -2  > 


(3.15) 


As  given  in  Appendix  H,  the  amplification  factor  |g|  for 
the  R-K  scheme  Eq.  (3.3a)  using  Eqs.  (3.12),  (3.14),  or 
(3.15)  is 


=  R  (CFL  ,  <J) )  +  kI(CFL,4>) 


(3.16) 
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where 


CFL  =  a  g  (3.17) 

<jj  =  wave  number  (linear  function  of  Ax) 

K  =  v^I 

and  R ( CFL ,  4> )  and  I(CFL,tj>)  are  the  real  and  imaginary  contri¬ 
butions  of  |G| ,  respectively,  and  are  somewhat  complicated 
functions  of  CFL  and  <J>*  The  R's  and  I's  resulting  from 
the  different  difference  representations  are  given  in 
Appendix  H. 

For  the  solution  to  remain  bounded  for  large  time, 
the  condition 


|G|  <  1  (3.18) 

must  prevail.  The  behavior  of  |g|  for  the  three  space- 
difference  approximations  used  is  shown  in  Figures  4a,  b, 
and  c  for  second-order  central,  first-order  backward,  and 
second-order  backward  differences,  respectively. 

Figure  4a  illustrates  the  desirable  feature  of  using 
second-order  central  differences  in  that  a  CFL  of  2/2 
(approximately  2.8)  can  be  used  (this  result  was  reported 
earlier  by  Jameson  et  al .  [47]);  whereas,  using  first-  and 
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|G| 


Figure  4.  Stability  Characteristics  for  Model  Problem 
Using  Four-Stage  Runge-Kutta. 
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second-order  backward  differences  require  using  CFL  numbers 
of  approximately  1.3  and  0.6,  respectively,  to  maintain 
stability  (Figure  4b  and  c) .  It  should  be  noted  that  the 
R-K  scheme  using  central  differences  is  stable  regardless 
of  the  sign  of  a  [45], 

Initially,  the  present  system  of  equations 
(Eq.  (2.41))  was  solved  using  central  space  differences 
with  a  CFL  =  2.8.  Smoothing  was  required  and  although 
reasonable  answers  were  obtained,  convergence  was  not  good 
and  the  amount  of  smoothing  needed  was  problem-dependent 
(smoothing  here  means  that  a  simple  weighted  averaging  was 
performed  after  each  time  cycle).  After  examining  the 
eigenvalues  of  the  NLj  and  matrices  in  Eq.  (3.9),  it 

was  found  that,  for  most  cases  considered,  eigenvalues  at 
each  grid  point  for  both  coefficient  matrices  were  positive. 
As  shown  by  Steger  et  al.  [52],  a  stable  scheme  for  the 
model  problem  (with  a  >  0)  can  be  constructed  using  back¬ 
ward  space  differences  which  has  better  dissipative  and 
dispersion  properties  than  that  of  a  centered  scheme. 

Figure  5  illustrates  the  convergence  obtained  for  the 
dummy  infinite  swept-wing  case  of  Cumpsty  and  Head  [16]  by 
solving  the  system  of  equations  (Eq.  (2.41))  using  second- 
order  central,  first-order  backward,  and  second-order 
backward  spatial  differences.  Here,  as  in  all  cases 
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Cycle 

Figure  5.  Convergence  Histories  for  the  Cumpsty  and  Head 
Infinite  Swept-Wing  Dummy  Test  Case. 
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computed,  convergence  is  measured  by  the  root-mean-square 
of  the  3H/0t  derivative,  that  is. 


NI,NJ  _  -> 

X  (I).  . 

i,j  1,1 

L  NI  •  NJ 


-1 1/2 


(3.19) 


For  the  solution  obtained  using  central  differences,  the 

following  simple  smoothing  scheme  was  used  for  the  solution 

—  T 

vector  U.  .  =  (0-,-.,  H,  A)  after  each  time-cycle: 

*vl  ]  ±  A. 


smoothed 


.  Hi+i.i+Hi-i 


.+u.  .,.+u.  .  ,+wU.  . 

,l~l, 1  +  1  ~i , 1-1  ~l ,  1 
uj  +  4 


(3.20) 


with  a)  =  8. 

From  Figure  5,  it  can  be  seen  that  convergence  for 
the  backward  schemes  using  no  smoothing  is  much  better 
than  that  obtained  using  central  differences  with  smoothing. 
The  convergence  obtained  using  the  backward  schemes  with 
no  smoothing  is  essentially  limited  by  the  machine 
truncation  error  (in  this  case,  approximately  14  signifi¬ 
cant  digits).  Also  shown  in  Figure  5  is  the  convergence 
history  of  the  first-order  backward  scheme  with  smoothing, 
which  illustrates  that  convergence  is  limited  when 
smoothing  is  applied. 

As  mentioned  earlier,  Jameson  et  al .  [47]  showed 
that  convergence  could  be  accelerated  by  using  local 
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time-steps  {recall  transient  results  were  not  considered 
important  here) .  Figure  5  also  illustrates  how  con¬ 
vergence  is  accelerated  for  the  present  system  of  equations 
by  using  spatially  variable  time-steps  as  opposed  to  using 
a  constant  time-step  (i.e.,  maximum  allowable  over  the 
field).  Therefore,  almost  all  of  the  solutions  presented 
herein  were  computed  using  first-order  backward  differ¬ 
encing  with  spatially  variable  time-stepping  and  no 
artificial  smoothing;  however,  as  discussed  in  the  next 
section,  central  differences  were  required  for  one  test 
case,  but  only  for  the  x^  spatial  derivative. 

3.3.  Boundary  and  Initial  Conditions 

Depicted  in  Figure  6  is  a  general  computational 
mesh  as  used  in  the  present  method.  The  boundary  con¬ 
ditions  used  to  compute  all  solutions  presented  here  were 
to  fix  or  "clamp"  the  three  dependent  variables  (0^,  H, 
and  A)  along  the  initial  start  line  (say,  along  the  leading 
edge  of  a  wing)  and  let  the  conditions  at  all  other 
boundaries  "float."  That  is,  the  dependent  variables  along 
all  boundaries  except  the  initial  start  line  were  computed 
in  the  same  fashion  as  those  in  the  interior  of  the  mesh 
using  appropriately  modified  spatial  differencing  near 
each  boundary.  This  essentially  amounts  to  extrapolating 
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□  -  Conditions  "Float"  }  Interior  Points 
X  -  Conditions  "Clamped"  )  Boundary 
O  -  Conditions  "Float'  i  Points 


Figure  6.  General  Computational  Domain- 
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the  conditions  at  the  boundaries  from  the  interior  of  the 
computational  domain. 

As  previously  mentioned,  all  eigenvalues  of  both 
the  and  matrices  in  Eq.  (3.9)  for  most  flow  cases 
are  positive  at  each  mesh  point.  Therefore,  specification 
of  boundary  conditions  along  the  initial  start  line  (i  =  1, 
j  =  1  -+■  NJ)  and  extrapolation  along  the  "outflow"  boundary 
(i  =  NI,  j  =  1  ■*  NJ)  is  compatible  with  the  sign  of  the 
characteristics  (eigenvalues).  However,  the  signs  of  the 
eigenvalues  of  the  matrix  were  mixed  at  some  mesh 

points  for  one  test  case  (to  be  presented)  and  central 
differences  were  used  to  approximate  the  D^fu1^)  deriva¬ 
tive  while  using  a  CFL  of  1.3  (see  Chapter  V,  Section  5.7). 
For  this  test  case,  extrapolation  of  boundary  conditions 
along  the  (i  =  1  -►  NI,  j  =  1)  line  and  ( i  =  1  ■+  NI,  j  =  N  J ) 
line  is  not  correct  if  one  adheres  strictly  to  the  infor¬ 
mation  obtained  from  the  sign  of  the  eigenvalues.  However, 
based  upon  the  results  to  be  presented,  it  seems  that  this 
erroneous  treatment  of  conditions  along  the  upper  and  lower 
boundaries  does  not  significantly  affect  the  outcome  of 
the  computations,  at  least  for  the  cases  considered. 

Solutions  generated  by  the  present  method  were 
found  to  be  insensitive  to  initial  conditions.  This  was 
determined  by  comparing  steady-state  solutions  obtained 
using  identical  edge  conditions  but  different  initial 
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conditions.  Therefore,  the  initial  conditions  used  for 
all  computations  presented  herein  were  to  set  the  values 
of  the  dependent  variables  at  each  mesh  point  to  those 
given  along  the  initial  start  line  (see  Figure  6). 


56 


AEDC-TR-83-37 


Chapter  IV 

COMPUTATION  OF  SURFACE  METRICS 

The  primary  objective  of  the  present  work  was  to 
develop  a  three-dimensional,  time-dependent,  turbulent 
integral  boundary-layer  computational  capability  that  can 
be  used  for  compressible  adiabatic  flow  in  nonorthogonal 

curvilinear  coordinates.  However,  before  any  boundary- 
layer  computations  can  take  place,  the  surface  metrics, 
i.e.,  h^  and  h 2  {scale  factors),  and  ,  K 2  an^  K21 

(curvature  terms)  must  be  known.  Most  of  the  cases  which 
have  been  computed  thus  far  were  such  that  specification 
of  the  surface  metrics  was  trivial  (i.e.,  =  h2  =  1  and 

=  K12  =  K21  =  0).  Therefore,  the  effects  of  these 
parameters  on  the  boundary-layer  calculations  in  a  general 
sense  are  not  included  in  the  present  study.  However,  two 
test  cases  have  been  calculated  which  use  surface  metrics 
other  than  those  for  the  trivial  case:  (1)  an  axisymmetric 
flow  where  the  metrics  are  determined  analytically,  and 
(2)  a  finite  swept  wing  case  using  a  skewed  mesh  where  the 
surface  metrics  are  computed  from  the  given  Cartesian 
coordinates . 

As  shown  in  Appendix  A,  all  metrics  depend 
explicitly  upon  derivatives  of  the  given  Cartesian 
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coordinates  with  respect  to  the  chosen  curvilinear  coordi¬ 
nates,  e.g.,  9x^/9x2-  For  the  axisymmetric  case,  these 
derivatives  are  evaluated  analytically;  whereas,  for  the 
finite  swept  wing  case,  the  surface  derivatives  were 
evaluated  using  simple  central  differences  except  near  the 
boundaries  where  extrapolation  was  used;  for  example, 

(h, }  =  (h, )  .  This  is  a  rather  crude  approxi- 

mation  compared  to,  for  example,  the  method  described  by 
Smith  and  Gaffney  [53]  who  approximate  the  metric  tensor 
using  bicubic  splines.  However,  computed  results  indicate 
that  the  approximations  used  here  were  adequate  for  the 
case  considered,  which  was  a  flat  surface. 
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Chapter  V 

RESULTS  OF  COMPUTATIONS 

Computations  using  the  present  method  are  compared 
with  the  results  of  seven  experimental/analytical  test 
cases  in  steady,  turbulent  flow.  These  data  sets  include: 
(1)  the  two-dimensional  (planar)  supercritical  RAE  airfoil 
flow  of  Cook,  McDonald,  and  Firmin  [54];  (2)  the  two- 
dimensional  (axisymmetric )  experiment  reported  by  Winter, 
Rotta,  and  Smith  [55]  concerning  the  flow  about  a  body-of- 
revolution;  the  infinite  swept-wing  cases  of  (3)  Cumpsty 
and  Head  [16],  (4)  van  den  Berg  and  Elsenaar  [56],  and 
(5)  Bradshaw  and  Terrell  [57]*  and  the  fully  three- 
dimensional  cases  of  (6)  East  and  Hoxey  [58],  and 
(7)  Humphreys  [25].  The  first  two  cases  were  run  to 
investigate  whether  the  results  of  the  present  three- 
dimensional,  time-dependent  method  duplicate  the  results 
of  its  steady,  two-dimensional  predecessor  [24];  whereas, 
the  remaining  cases  test  the  three-dimensional  compu¬ 
tational  capabilities  of  the  code.  In  addition,  results 
using  the  present  method  are  compared  with  results  of  other 
calculations  when  available. 
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5.1.  RAE  Airf oi 1 --Two-Dimensional  Planar 

The  experiment  of  Cook  et  al.  [54]  involved  the 
measurement  of  turbulent  boundary-layer  quantities  over  a 
planar,  supercritical  airfoil  in  adiabatic,  high  Reynolds 
number,  transonic  flow  at  various  angles  of  attack.  The 
particular  flow  conditions  at  which  comparisons  will  be 
made  are 


M  =  0.600 

CO 

Re  =  6 . 3xl06 
00,  c 

and 

2.57  deg  angle  of  attack 


which  is  denoted  as  Case  3  in  [54]. 

Comparisons  between  experiment  and  results  computed 
by  the  present  method  are  shown  in  Figure  7  for  displace¬ 
ment  and  momentum  thicknesses,  shape  factor,  and  skin 
friction.  Also  given  in  Figure  7  are  the  results  of  the 
calculation  method  described  in  [24].  Comparing  the 
results  of  the  present  method  with  those  of  [24]  reveals 
that  the  curves  are  virtually  indistinguishable,  except 
for  shape  factor  near  the  airfoil  trailing  edge,  where  the 
present  results  are  slightly  lower.  In  addition,  the 
agreement  between  measured  and  calculated  results  is  con¬ 
sidered  good. 


60 


[5IS1 


AEDC-TR-8  3-3  7 

5.2.  Waisted  Body  of  Revolution — 
Two-Dimensional  Axi symmetric 


The  objective  of  the  investigation  reported  by 
Winter  et  al.  [55]  was  to  produce  an  axisymmetric,  con¬ 
verging  flow  with  an  adverse  pressure  gradient  to  determine 
the  effects  of  Mach  number,  pressure  gradient,  and  stream¬ 
line  convergence  and  divergence  on  adiabatic,  turbulent 
boundary-layer  development.  Comparisons  between  results 
computed  by  the  present  method  and  by  the  method  in  [24] 
with  experiment  are  made  at  the  following  freestream 
conditions : 


M  =  0.597 

CO 

Re  _  =  9 . 98xl06 

00  |  Jj 


The  present  computations  were  made  using  the  following 
metric  coefficients: 


h,  =  1 


h„  = 


r T  =  local 
w 


K-,  =  0 


K„  =  - 


1  drw 


rw  dx 


1L 


K 


12 


K 


21 


=  0 

:_K, 


body  radius 


1/2 
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Figure  8a  illustrates  measured  and  computed  distri¬ 
butions  of  skin  friction,  displacement  and  momentum 
thickness,  and  shape  factor,  and  the  agreement  is  con¬ 
sidered  good.  Differences  between  the  present  computations 
and  those  of  [24]  are  seen  to  be  rather  small.  In  addition. 
Figure  8b  gives  comparisons  between  measured  and  computed 
boundary-layer  velocity  profiles  at  five  axial  locations 
down  the  body;  the  agreement  is  considered  excellent. 

The  results  presented  in  Figures  7  and  8  for  these 
special  cases  (i.e.,  two-dimensional-planar  and 
axi symmetric )  lend  credence  to  the  present  method  in  that 
results  generated  by  a  previous  method  were  essentially 
duplicated.  The  remaining  cases  were  chosen  such  that  the 
three-dimensional  capabilities  of  the  present  method  could 
be  investigated. 

5.3.  Cumpsty  and  Head--Inf ini te  Swept  Wing 

The  Cumpsty  and  Head  [16]  "dummy"  test  case  is 
that  of  an  infinite  swept  wing  in  an  incompressible  flow 
with  a  constant  linear  gradient  of  chordwise  velocity  and 
a  constant  spanwise  velocity,  given  by 
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%  =  0.597 
Reoo,  L  =  9-98  x  10,5 


a  o  Measured  f  55] 


a.  Skin  Friction,  Momentum  and  Displacement 
Thicknesses,  and  Shape  Factor 


us'q 


b.  Velocity  Profiles 


Figure  8.  Computed  and  Measured  Boundary-Layer  Quantities 
for  the  Waisted  Body  of  Revolution  [55]. 
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(  q  cos  a  ,  x_  <  0 
!  o  '  c 


cos  aQ  (1  -  kxc)  f  xc  > 


u 

e 


s 


q  sin  a 
^oo  O 


where  ue  and  ue  are  chordwise  and  spanwise  edge  veloci- 
c  s 

ties,  respectively,  aQ  is  the  sweep  angle,  and  k  is  the 

velocity  gradient.  The  cases  considered  here  are  for 

a  =  35  deg,  k  =  0.250,  and  k  =  0.267.  All  calculations 
o 

were  begun  at  xc  —  0 . 

Comparisons  of  results  between  the  present  method 
and  the  calculations  of  Cumpsty  and  Head  [16J  for  momentum 
thickness,  shape  factor,  wall  streamline  angle,  and  chord- 
wise  skin  friction  are  shown  in  Figure  9  for  the  case  of 
k  =  0.267.  The  agreement  between  the  two  calculation 
methods  is  considered  good.  It  should  be  noted  that  the 
computations  of  Cumpsty  and  Head  shown  in  Figure  9  are 
those  resulting  from  using  all  the  terms  in  their  equations 
with  overall  iteration  (curves  labeled  "7(a)"  in  [16]). 
Figure  10  compares  the  present  calculations  of  chordwise 
skin  friction  with  those  of  Cumpsty  and  Head  [16], 

Cebeci  [59],  and  Bradshaw  [60]  for  the  k  =  0.250  case. 
Although  the  present  computations  are  slightly  lower  than 
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a.  Momentum  Thickness 


1  2  h 

i  o  r  i  i  i  i  i  i  i  i _ i  i  i  i  i  i 

0  0.2  0.4  0.6  0.8  1.1  1.3  1.5 

xc,  ft 


b.  Shape  Factor 


d.  Chordwise  Skin  Friction 


Figure  9.  Calculated  Results  for  the  Dummy  Test  Case  of 
Cumpsty  and  Head  [16], 
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Range  of  Values  Computed  by  Cebeci  [59] 
Bradshaw  [60]  and  Cumpsty  and  Head  [16] 


Figure  10.  Calculated  Skin  Friction  Distribution  for  the 
Dummy  Test  Case  of  Cumpsty  and  Head  [16]. 
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those  of  the  other  investigators,  correct  trends  are 
indicated. 


5.4.  van  den  Berg  and  Elsenaar--Inf inite  Swept  Wing 

The  low  speed  flow  of  van  den  Berg  and  Elsenaar 
[56]  involved  probing  the  three-dimensional  boundary  layer 
on  a  flat  surface  swept  at  35  deg  with  an  external  pres¬ 
sure  distribution  induced  by  an  appropriately  shaped  body 
such  that  infinite  swept-wing  conditions  were  approxi¬ 
mately  simulated  (this  experiment  was  performed  specifi¬ 
cally  for  comparison  to  computational  methods).  Figure  11 
illustrates  the  measured  flow  angle  and  surface  pressure 
distributions.  Also  shown  in  Figure  11  are  the  analytic 
representations  of  measured  flow  angles  and  pressures.  As 
suggested  in  [56],  the  analytical  representations  of  the 
external  flow  conditions  were  used  as  input  in  the  present 
computations.  It  is  seen  that  measured  pressures  are 
represented  well  over  the  entire  axial  distance,  whereas 
the  expression  for  flow  angle  falls  below  those  measured 
for  >  1.05  m.  Comparisons  between  measured  and  computed 
boundary-layer  quantities  are  shown  in  Figure  12.  Agree¬ 
ment  between  the  computations  and  measurements  upstream  of 
x^  =  1  meter  is  considered  good;  whereas,  past  this  point, 
considerable  discrepancies  exist,  particularly  with  the 
streamwise  integral  thicknesses.  Similar  results  were 
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20 


b.  Surface  Pressures 

Figure  11.  Inputs  for  van  den  Berg  and  Elsenaar  [56] 
Infinite  Swept-Wing  Case. 
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xr  m 


c.  Wall  Streamline  Angle 

Figure  12.  Computed  and  Measured  Boundary-Layer  Quantities 
for  the  van  den  Berg  and  Elsenaar  [56] 

Infinite  Swept-Wing  Case. 
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obtained  by  Cebeci  and  Chang  [61]  using  a  finite-difference 
method.  However,  as  shown  in  Figure  11a,  the  measured 
flow  angles  downstream  of  =  1  meter  deviated  from  those 
which  were  used  as  input  for  the  computations  (the 
analytical  representation  was  that  for  an  infinite  swept- 
wing  [56]).  Illustrated  in  Figure  13  is  a  comparison 
between  computed  wall  streamlines  and  a  surface  oil  flow 
photograph  which  indicates  reasonable  qualitative  agree¬ 
ment  between  the  experiment  and  the  computations. 

An  unsuccessful  attempt  was  made  to  compute  this 
flow  using  the  measured  flow  angles  as  given  by  the  symbols 
in  Figure  11a.  According  to  the  experiment  [56],  the  flow 
was  separated  near  the  trailing  edge  and  it  is  suspected 
that  this  phenomenon  is  responsible  for  an  instability  in 
the  code  which  leads  to  program  failure.  (In  fact,  Delery 
and  Formery  [62]  use  this  experiment  to  test  their  inverse 
method  which  also  indicates  that  the  flow  was  separated.) 

5.5.  Bradshaw  and  Terrell — Infinite  Swept  Wing 

The  experiment  of  Bradshaw  and  Terrell  [57]  was 
set  up  to  test  the  outer-layer  assumptions  in  extending 
Bradshaw's  et  al.,  calculation  method  [63]  from  two  to 
three  dimensions.  This  was  an  infinite  swept  wing  in  a 
low-speed  flow  where  boundary-layer  measurements  were 
obtained  over  the  rear,  flat  portion  of  the  wing.  The 
axial  pressure  gradient  was  nominally  zero  with  a  decaying 
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crossflow.  Comparisons  of  computed  and  measured  boundary- 
layer  quantities  are  given  in  Figure  14.  Agreement  among 
measured  and  calculated  skin  friction,  cross-flow  angle 
distribution,  and  wall  streamline  angle  is  not  as  good  as 
the  computations  of  Cebeci  [59]  and  Bradshaw  [60]  (which 
are  differential  methods).  For  this  flow,  East  [10] 
reported  similar  findings  regarding  the  relative  perform¬ 
ance  of  integral  and  differential  methods,  and  blamed  the 
poor  performance  of  integral  methods  on  the  use  of 
Johnston’s  [29]  cross-flow  profile  which  is  not  repre¬ 
sentative  of  a  decaying  three-dimensional  flow,  particu¬ 
larly  in  the  wall  region.  East  [10]  points  out  that  this 
flow  demonstrates  how  differential  methods  are  generally 
more  flexible  than  integral  methods  when  applied  to  a 
variety  of  different  types  of  flows.  However,  as  shown  in 
Figure  14c,  the  present  computations  of  streamwise  velocity 
agree  very  well  with  experiment. 

5.6.  East  and  Hoxey — Fully  Three-Dimensional 

The  experiment  of  East  and  Hoxey  [58]  consisted  of 
an  obstruction  placed  in  a  thick  two-dimensional  boundary 
layer  on  the  floor  of  a  low-speed  wind  tunnel.  Pressure 
gradients  caused  by  the  obstruction  induced  three- 
dimensionality  and  separation.  Boundary- layer  quantities 
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Figure  14.  Computed  and  Measured  Boundary-Layer  Quantities 
for  the  Bradshaw  and  Terrell  [57] 

Infinite  Swept-Wing  Case. 
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were  measured  upstream  of  the  obstruction  in  a  region  off 
the  plane  of  symmetry.  Measured  wall  pressures  and  edge 
flow  angles  were  used  as  inputs  to  the  computations. 

Figure  15  gives  comparisons  of  measured  and  computed 
boundary-layer  parameters  6  inches  from  the  plane  of 
symmetry.  Also  shown  are  the  results  of  Myring's  integral 
method  [13].  Good  qualitative  trends  are  obtained  with 
the  present  method,  although  values  of  momentum  thickness, 
shape  factor,  and  wall  streamline  angle  are  generally 
underpredicted  for  x-^  <  25  inches.  Agreement  between 
measured  and  calculated  streamwise  velocity  profiles  is 
considered  satisfactory. 

Illustrated  in  Figure  16  are  computed  streamline 
patterns  at  the  wall  (those  at  the  edge  are  input)  along 
with  those  as  deduced  from  the  flow  field  measurements. 
Reasonable  qualitative  trends  are  indicated  by  the  compu¬ 
tations  except  near  the  separation  line  (indicated  in 
Figure  16b),  which  is  not  surprising  because  the  present 
computations  did  not  predict  separation. 

5.7.  1978  Stockholm  Test  Case--Finite  Swept  Wing 

The  1978  Stockholm  Test  Case  [25]  was  based  upon 

the  flow  about  a  swept  wing  of  modern  configuration  in 

6 

high  subsonic  flow  (M^  =  0.5,  Re^  =  7x10  /unit  length). 
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o  Measured  [58] 
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Figure  15.  Computed  and  Measured  Boundary-Layer  Quantities 
for  the  Experiment  of  East  and  Hoxey  [58]. 
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Figure  16.  Computed  and  Measured  Wall  Streamlines  for  the 
East  and  Hoxey  Test  Case  [58], 
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Effects  of  three-dimensionality  and  compressibility  were 
small  but  non-negl igible .  Because  the  test  case  was 
designed  to  test  three-dimensional  boundary- layer  calcu¬ 
lation  methods,  the  inviscid  velocity  distribution  was 
provided  by  a  higher  order  panel  method-  The  only  experi¬ 
mental  data  provided  for  comparison  were  oil  flow  photo¬ 
graphs  of  the  wing  surface.  Boundary  values  at  the  initial 
start  line  are  given  in  [25].  Figure  17  gives  the  results 
at  span  stations  2,  4,  and  6  as  computed  by  the  present 
method  along  with  the  range  of  all  eight  calculation 
methods  which  were  compared  in  [25]  for  the  case  of  zero 
angle  of  attack.  Favorable  agreement  is  seen  to  exist 
when  comparing  momentum  thickness,  shape  factor,  and  skin 
friction;  whereas,  the  present  computations  for  Bw  indicate 
a  more  rapid  increase  for  x/c  >  0.4  at  all  span  stations 
than  do  the  other  calculations.  Figure  18  gives  the  wall 
streamline  patterns  as  computed  by  the  present  method  and 
as  measured  using  oil  flow  visualization;  good  qualitative 
agreement  is  seen  to  exist  between  the  computations  and 
the  measurements. 

Similar  to  the  van  den  Berg  and  Elsenaar  infinite 
swept-wing  case  [56]  discussed  in  Section  5.4,  another 
program  failure  was  encountered  when  an  attempt  was  made 
to  compute  the  8-deg  angle-of-attack  case  reported  in  [25]. 
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Figure  17.  Calculated  Boundary-Layer  Parameters  for  the 
1978  Stockholm  Test  Case  [25]  for 
Zero  Angle  of  Attack. 
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Although  the  general  consensus  from  the  results  reported 
in  [25]  was  that  the  flow  was  not  separated  [also  implied 
by  surface  flow  visualization),  it  is  again  suspected  that 
separation  is  responsible  for  the  observed  program  failure. 

For  this  particular  test  case,  the  boundary  con¬ 
ditions  were  not  treated  correctly  according  to  the  sign 
of  the  characteristics  along  both  "side"  boundaries;  that 
is,  (i  =  1  +  NI,  j=l)  and  (i  =  1  -+  NI,  j  -  NJ )  (refer  to 
Figure  6,  page  54).  The  results  presented  in  Figures  17 
and  18  indicate  that  this  erroneous  treatment  of  boundary 
conditions  seems  to  have  had  little  effect  on  the  compu¬ 
tations  for  the  zero  angle-of-attack  case.  However,  as 
noted  above,  the  numerical  procedure  became  unstable  when 
attempting  to  compute  the  8-deg  angle-of-attack  case;  what 
effect  (if  any)  the  boundary  conditions  along  the  two 
"side"  boundaries  had  on  the  computations  for  this  test 
case  is  not  known. 
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Chapter  VI 

SUMMARY,  CONCLUSIONS,  AND  RECOMMENDATIONS 

A  method  for  computing  three-dimensional,  time- 
dependent,  compressible,  turbulent  boundary  layers  in 
nonorthogonal  curvilinear  coordinates  has  been  presented. 
The  method  solves  the  time-dependent  momentum  and  mean- 
flow  kinetic  energy  boundary-layer  integral  equations 
which  provides  the  viscous  portion  of  a  viscous/inviscid 
interaction  approach  where  identical  surface  grids  for 
both  the  viscous  and  inviscid  calculation  methods  can  be 
used.  A  four-stage  Runge-Kutta  time-stepping  scheme  was 
used  to  numerically  solve  the  system  of  equations  using 
local  time  steps  to  accelerate  convergence.  Several  space 
difference  approximations  were  employed  and  it  was  found 
that  a  backward  scheme  gave  the  best  results  using  no 
artificial  smoothing.  Calculated  results  using  the  present 
method  were  compared  to  experimental  data  and  to  results 
of  other  calculation  methods;  satisfactory  agreement  was 
obtained . 

There  are  several  areas  where  improvements  could 
be  incorporated  which  would  in  turn  improve  the  overall 
outcome  and  performance  of  the  present  computational  pro¬ 
cedure.  Firstly,  incorporation  of  a  more  general  method 
of  computing  surface  metrics  is  needed  in  order  to  properly 
address  three-dimensional  geometries;  for  example,  the 
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procedures  of  Smith  and  Gaffney  (53)#  or  Craidon  [64]. 

The  use  of  a  more  general  cross-flow  velocity  profile 
representation  would  make  the  present  method  more  flexible 
with  respect  to  the  classes  of  flows  which  could  be 
addressed.  Although  other  models  are  available  (e.g., 
Mager's  cross-flow  model  [65]),  none  have  sufficient  flexi¬ 
bility  to  allow  a  three-dimensional  integral  method  to  be 
routinely  applied  to  any  flow  field  such  that  accurate 
results  can  generally  be  expected.  As  discussed  in 
Appendix  D,  a  more  thorough  treatment  of  the  "cross-flow” 
dissipation  is  needed.  However,  to  properly  handle  this 
term  in  the  equations  requires  an  accurate  cross-flow 
model  which  is  valid  for  a  wide  variety  of  three- 
dimensional  boundary-layer  flows. 

All  calculations  presented  herein  were  performed 
using  a  CRAY-1S  computer.  Depending  on  the  number  of  mesh 
points  and  type  of  spatial  differences  used,  solution 
times  ranged  from  15  to  300  CPU  seconds.  This  is  somewhat 
slow  considering  the  relative  simplicity  of  the  geometries 
involved.  However,  it  is  felt  that  performance  of  the 
numerical  procedure  could  be  improved  significantly  by  more 
efficient  coding  and  by  incorporating  a  different  numerical 
scheme.  For  example,  a  two-stage  R-K  scheme  operating  at 
a  CFL  of  0.9  was  incorporated  into  the  code  and  the 
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van  den  Berg  and  Elsenaar  case  [56]  was  recomputed.  The 
solution  was  practically  identical  to  that  using  the  four- 
stage  scheme  with  a  40%  reduction  in  CPU  time. 

It  is  apparent  that  the  major  deficiency  of  the 
present  method  is  its  inability  to  compute  some  flow 
fields:  in  the  present  study,  the  van  den  Berg  and 
Elsenaar  case  [56]  (infinite  swept  wing)  using  the  measured 
wall  pressures  and  edge  flow  angles,  and  the  1978  Stockholm 
test  case  [25]  (finite  swept  wing)  at  8  deg  angle  of 
attack.  This  seems  to  be  a  significant  hindrance  when 
considering  using  this  method  as  one  part  of  a  viscous- 
inviscid  interaction  approach.  On  the  other  hand,  knowing 
that  the  method  has  failed  for  a  particular  case  could 
possibly  be  interpreted  as  signaling  whether  or  not  sepa¬ 
ration  has  been  encountered,  although  no  information  is 
generated  pertaining  to  its  location.  As  discussed  by 
Cousteix  and  Houdeville  [66],  and  Delery  and  Formery  [62], 
it  is  apparent  that  an  inverse  formulation  of  the  boundary- 
layer  equations  is  required  if  separated  flows  are  to  be 
addressed,  even  for  a  time-dependent  computational  method. 
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APPENDIX  A 


SURFACE  METRIC  COEFFICIENTS 


Before  boundary-layer  computations  can  proceed, 
metrics  associated  with  surface  scale  factors  and  curvatures 
must  be  determined.  The  following  relationships  for  a 
general,  nonorthogonal ,  curvilinear  coordinate  system 
embedded  in  the  surface  of  a  body  (see  Figure  1,  page  8) 
are  listed  by  Cebeci  et  al,  117]  and  are  given  here  for 
completeness  in  slightly  different  form.  As  before,  the 
subscript  "i"  takes  on  values  of  1  and  2,  and  subscript  3. 
resulting  from  i  +  1  when  i  =  2  is  taken  to  be  subscript  1. 


(A.l) 


K, 


Ki,i+1 


1  r  1  fa  . 

h-^h2  sinA  JJijl  3x^ 


hi+1  cos* 


1 

sinX 


+  cosX 


(A. 2) 


(A. 3) 


(A. 4) 
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(A. 7) 


(It  should  be  pointed  out  that  h. ,  K . ,  and  K.  .  .  must 

11  1  /  1 ' 1 

be  computed  at  each  mesh  point  in  the  computation  domain.) 
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APPENDIX  B 


DEFINITION  OF  STREAMLINE  INTEGRAL  LENGTHS 


Using  the  sketch  below,  velocities  in  the  non- 
orthogonal  axis  system  are  related  to  those  in  a  streamline 
system  by  the  following  relations  [13,14]: 


u 


1 


u  sin£  -  u  cost, 
s _ _ _ n 

■ 


(B.l) 


u 


2 


u  sina  +  u  cosa 
s _ n _ 

sinX 


(B.2) 


At  the  boundary-layer  edge,  =  0.  Thus, 


u 


1 


—  sm£ 

U  — T - T- 

s  sinA 


( B ,  3 ) 


u 


2 


—  sina 
u  — t — r 
s  smX 


(B,  4) 


where  (u^,  U2)  are  velocities  resolved  in  the  nonorthogonal 

system,  and  (u  ,  u  )  are  those  in  the  streamline  system 
3  n 

(recall  overbars  denote  values  evaluated  at  the  boundary- 
layer  edge)  . 
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Xl-  U1 


Sketch  Showing  Relationship  Between  Streamline 
and  Nonorthogonal  Coordinates  (Planform  View) 


Also,  the  resultant  velocity  is  given  by 

2  2  2  2  2 
q  =  u.  +  u.  +  2u,  u~  cosA  =  u  +  u  (B.5) 

12  12  s  n 

and  at  the  boundary-layer  edge, 

—2  —2  —2  —  —  —2 
q  =  u,  +  u~  +  2u,u~  cosA  =  u  (B.6) 

^12  12  S 

Thus,  integral  lengths  in  the  streamline  coordinate  system 
can  be  defined  as  [13,14]: 


(B .  7 ) 

(B.8) 

(B.9) 
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P  q2  ©12 

p  q2  ©21 

P  q  022 
P  q3  El1 

p  q3  e12 
p  ?3  e21 

P  q3  E22 

q  A* 

^  u 

_  * 

q  A 

V 


/  Pun<us-Us)dx3 


‘  /  pusun  dx3 


/  pu  dx. 
0  n  3 


0 


;  Pus  ( -  u3  )  dx3 


0 


1  pun  (  uj  -  u2s  )  dx3 


-  /  pusun  dx3 


-  /  pu  dx_ 

0  n  j 


/  (“s-us)dx3 


/  dx 
0  n  3 


(B. 10) 

(B.ll) 

(B . 12) 

(B . 13) 

(B . 14 } 

(B . 15) 

(B ,  16) 

(B . 1 7 ) 

(B.18) 
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APPENDIX  C 

STREAMWISE  SHAPE  FACTOR  CORRELATIONS 

Except  for  the  correlation  derived  by  Donegan  [37] 

for  H0  ,  all  shape  factor  correlations  used  in  the  present 
P 

study  are  listed  in  [24] ,  [42] ,  and  [43] ,  and  are  given 

here  for  completeness: 

H  =  H  (  1  +0,113  M2  )  +  0,29  M2  (C.l) 

6  © 


(He*) 


Me=0 


=  1.48061  +  3.83781e 


— 2H 


+  tan_1 


107"H-1 


L  1.23 


-  (.  0.33  T’T~T  )  tanh1/2 
1  /  ♦  i 


/  _6W  7.HN1-45761  ] 

(  1.2874x10  b  )  (  10  )  J 


(C .  2 ) 


He*  = 


(h9*)m  =o  +  °*028  K 

_ e _ 

1  +  0.014  M2 


(C.  3 ) 


0 


0.92  M 


11  =  1  - 

Gu  7.09  +  M 


j  tanh [1 .49 (H  -  0.9)]  (C.4) 


He  =  M*  (0,185  H  +  0.150) 
P  e 


(C.5) 
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c  DS 
f  u 


0 . 01167e  °'038h3  +  (  9.0xl0“8)e 


-8\ J 


where 


+  0.0115  +  ACFD/1000 


/(l  +  ( 


ACFD  =  mRe 


G 


11 


m  =  650H  -  743 


n  =  -  1. 59H  +  1. 45  J 


ACFD  =  menRe0n 


m  =  3.25  e 


0.045H2 


n  =  H/10000  -  0.0017 


..  603H 

1.025  Me1,4  )  (C .  6} 

H  £  1.6 

H  >  1.6 
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APPENDIX  D 

DISSIPATION  INTEGRALS 

The  resolution  of  the  dissipation  integrals 
(appearing  as  the  last  term  in  the  mean-flow  kinetic 
energy  integral  equation,  Eq.  (2.20))  from  quantities  in 
nonorthogonal  coordinates  to  those  in  the  streamline 
coordinate  system  is  given  in  this  Appendix.  These 
manipulations  are  performed  such  that  correlations  derived 
for  two-dimensional  turbulent  boundary  layers  can  be 
utilized  for  the  general  three-dimensional  case. 

From  Eq.  (2.20),  the  dissipation  integrals  are 
written  as 


D  +  D 
U1  U2 


/ 

0 


9t 


X1 


3x. 


+  u. 


(D.l) 


Integrating  by  parts  and  using  the  boundary  conditions 
that 


It  J  =  I  u.  J  -  [  u  ]  =0  (i  =  1  or  2)  (D. 2) 

i  °°  0  z  0 


Equation  (D.l)  becomes 


9u. 


9u, 


Du  +  Du  “ _ 3  f  I  Tx  'ot;  +  tx  - — —  I  dxQ 

1  2  p  q  0  \  X1  9x3  2  3x3  I  3 


(D.3) 
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Using  the  relations  between  velocities  in  nonorthogonal 
and  streamline  coordinates  (see  Eqs.  (B.l)  and  (B.2)  in 
Appendix  B) ,  Eq.  (D.3)  becomes 


D  +  D  -1 

ui  u2  -^3 — 7  , 

p  q  sinA 


00  3u 


■sin? 

/ 

L 

0 

00 

3u 

;  Tx 
X1 

I 

3x 

00 

3u 

/  Tx 
X2 

3x 

00 

3u 

1  T*0 

1 

3x. 

n 


[x,  w:  ^3 


dx. 


n 


dx. 


dx. 


(D.4) 


Substituting  the  outer  portion  (Eq.  (2.31b))  of 
Johnston's  cross-flow  profile  [29]  into  the  above  yields 


D  +  D 
U1  u2 


p  q3  sinA 


(sin£+Acos£)  I  x 
0 


3us  , 

X1  *3  3 


<»  3u 


+  (s inct’-Acosa)  /  t 


0  x2  3x3 


] 


(D.5) 


The  total  shear  stress  in  the  and  Xj  directions  can  be 
resolved  into  streamwise  components  as 
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Txx  =  iIH3T  <  Ts  sin«  -  Tn  cos5  1 


(D.  6a) 


X  =  — : - r-  (  T  Sind  +  T  COSOt  ) 

x2  sinA  s  n 


(D. 6b) 


where  xg  is  the  streamwise  component  of  the  total  shear 
stress  and  x^  is  normal  to  it.  Substituting  Eq.'  (D.6) 
into  (D.5)  gives 


D  +  D 


,  C  *  DS  /  Dn 

1  f  U  L  A  t  u 

.  2  5  ri  t2  _s 

sin  X  \  D 

'  u 


where 


(D.7) 


tl  = 


t2 


2  .  2 

sin  £  +  sm  a+A(sin£  cos£  ~  sina  cosa) 

.  .  2  2 
sin£  cos^  +  sincx  cosa-A(cos  £  +  cos  a) 


(D.8a) 


{D. 8b) 


Ii_  JV2 1  dx 
0  %  3 


<D.8c) 


.  x  m 
0;  tanB  v  ~&r dX3 

w 


(D.8d) 


and  the  relation 


x  =  x  tang 
ns 


(D.9) 


has  been  used  in  Eq.  (D.8d) 
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As  discussed  in  Section  2.2.4,  Eq.  (D.7)  is  evalu¬ 
ated  in  the  present  study  using  a  correlation  for 
c  D®/2  [42,43]  and  assuming  that 


D 

D 


n 

u_ 

s 

u 


<  < 


1 


(D. 10} 


In  principle,  the  "crossflow"  dissipation  (Eq.  (D.8d)) 

could  be  correlated  in  an  analogous  manner  as  was  done  for 

"streamwise”  dissipation  if  tang  could  be  evaluated 

analytically  for  the  general  case.  For  example,  the  total 

shear  stress  components  in  streamline  coordinates  are 
1 

given  as 


T 


s 


pu'u- 


(D.lla) 


T 

n 


pu  1  u ' 
K  n  3 


( D. lib) 


■*"The  overbars  and  primes  denote  time-averaged  and 
fluctuating  quantities,  respectively;  e.g., 


I  u.  _ 

s  3 


1 

2T 


t-T 


t+T 

/  (u^)dt 


where  the  averaging  period  T  is  much  larger  than  the  time 
scale  of  the  fluctuations  (see  [31]  and  [44]). 
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The  Reynolds  stresses  can  be  modeled  using  the  eddy 
viscosity  concept  (e.g.,  see  [59])  as 


-  pu '  u  ’ 

S  X, 


(D.12a) 


9u 


"  punV 


n 


9x. 


(D.12b) 


where  p. 


molecular 


and  p  are  the  eddy  viscosities,  and  p  is  the 
r2 

viscosity.  Thus,  Eqs .  (D.lla,b)  become 


T 

S 


(  y  + 

1 


T 


n 


(  ij  +  if 
U2 


Therefore,  tanB  could  be  evaluated  using 

T 

tanB  =  — 

T  _ 


(D  .13a) 


(D. 13b) 


(D.14) 


by  using  Johnston's  cross-flow  profile  [29]  for  un  and  by 

making  appropriate  assumptions  concerning  the  evaluation 

of  the  eddy  viscosities  p  and  p  (e.g.,  see  [59]). 

fcl  t2 
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However,  similar  to  the  "streamwise"  dissipation 
correlation  [42,43]  ,  which  is  based  upon  a  very  general 
description  of  the  streamwise  velocity  profile  [28]  and  an 
eddy  viscosity  formulation  (yt  )  derived  from  two- 
dimensional  boundary-layer  analysis  [44],  a  correlation  of 
the  "crossflow"  dissipation  would,  in  turn,  be  based  upon 
the  cross-flow  velocity  profile  plus  an  additional  eddy 
viscosity  formulation  (pt  ) .  Therefore,  in  lieu  of  the 
lack  of  generalities  of  both  the  cross-flow  profile  repre¬ 


sentation  and  the  eddy  viscosity  ,  the  contrit)ut:i-on 

was  simply  neglected  in  comparison  with  D® .  Obviously, 
more  study  is  needed  concerning  the  evaluation  of  . 
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APPENDIX  E 


RELATIONS  BETWEEN  THE  STREAMWISE  INTEGRAL  LENGTHS 


For  a  given  cross-flow  velocity  profile,  certain 
relations  exist  between  streamwise  integral  lengths  which 
enable  the  13  unknowns  to  be  expressed  in  terms  of  five. 
These  relations  are  given  here  for  the  Johnston  cross- 
flow  profile  [29]  only. 

The  13  unknowns  in  streamline  coordinates  are 


A1 '  A2'  Qll' 


012'  °21'  022 '  Ell'  E12 '  E21 '  E22'  Au' 


and  0^  (see  Appendix  B  for  these  definitions).  As  an 
example  of  how  these  integral  lengths  are  interrelated, 
consider  Eq.  (B.12)  : 


- 2 

p  q 


/  pun  dx3 


(E.l) 


Thus,  using  only  the  outer  portion  of  Johnston's  cross- 
flow  model  [29,14] f 


P  q2  © 


22 


A2  /  (  pu2-2pusus+pu2  )  dx3 


2  —  —  —  —2—  — 

A  /  [  us  (  p  us-pus  )  -us(p-p)-pug  (  us-us) ]dx3 


A2  t  p  uj  A*-p  u2  ep-p  u2  0U  ) 
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or, 

®22  ~  A  ^  ^11  +  )  (E.2) 

The  remaining  relationships  can  be  derived  in  an  analogous 
manner,  and  are  summarized  as  follows: 


A*  =  -  A  (  A*  -  ep  )  (E.3) 


012 

=  A  (  A*  -  90  -  eu  ) 

(E.4) 

021 

=  -  A  0U 

(E.5) 

°22 

-  <  0n  +  «P  -  > 

(E.6) 

E21 

=  ft2  (  En  -  2en  ) 

(E.7) 

E12 

=  A  (  Aj  +  0U  -  En  -  6b  ) 

{E .  8) 

E22 

A  (  A1  *  +  EH  "  3®11  1 

(E  ■  9 ) 

a; 

=  -  A  A* 
u 

(E.10) 

It  is  computationally  more  efficient  to  use 

the  shape 

factors  as  defined  in  Eq.  (2.27)  and  rewrite 

(E.10)  as 

Eqs .  (E.3)  to 

021  “  A  011 

(E. 11) 

a2*  =  e2i  (H  -  He  > 

p 

(E.12) 
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12 

021  - 

A  * 

A2 

(E. 13) 

22  = 

~  A  012 

(E, 14 ) 

12 

G12  + 

G21  (H0*  -  2  > 

(E  .15) 

‘21  = 

-  022 

A  E12 

( E . 16 ) 

'22  = 

-  A  (  E21 

.  ~  °22  > 

{E .17) 

A*  = 

V 

-  A  A* 

(E .17) 

no 
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APPENDIX  F 

RELATIONS  BETWEEN  INTEGRAL  LENGTHS  IN  STREAMLINE 
AND  NONORTHOGONAL  COORDINATE  SYSTEMS 

The  derivation  of  the  relationships  between  the 
integral  thicknesses  in  the  two  coordinate  systems  is 
algebraically  tedious  but  straightforward  (as  cited 
earlier,  Smith  [14]  gives  as  an  example  the  derivation  of 

911  in  terrns  of  Qn'  012'  etc-)-  The  complete  list  of  the 
12  relationships  is  as  follows  (recall  that  upper-case 
Greek  letters  denote  integral  lengths  using  velocities 
expressed  in  the  streamline  coordinate  system;  whereas, 
lower-case  Greek  letters  represent  those  resolved  in  the 
nonorthogonal  system) : 


sinX  (41  Sln«  -fl2  cos5) 

(F.l) 

•i  - 

\  it  * 

— : — t-  (A.  sina+A_  cosa) 
sinA  1  2 

(F.2) 

en  = 

— [0X1  sin2^  -  (012  +  ©21)sinC 

cos? 

sin  A 

2 

+  022  cos  £] 

(F.3) 

S12  - 

— [Qti  sina  sin£  +  012  sinf;  cosa 

sm  a 

-  02^  cos£  sina-022  cosa 

cosEJ 

(F.4) 

II 

CD 

—  ~*~2'  [0,  i  sina  sin£  +  0^  sini?  cosa 

sin  a 

-  0^2  cos£  sina  -  ©22  cosa 

cos£] 

(F .  5) 
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22 


1  2 

- Y~  [0-q  sin  a  +  (012  +  ©2i)  cosa  sina 


sin  A 


+  ©22  cos  a] 


'11 


I  ”3  o 

— xz  ^Eii  sin  ?  “  (3E12  +  2A*)  sirt  ?  cos? 


sin3*  1  1J- 


+  3E21  sin?  cos2?-E22  cos3?] 

1  O  JL 

[E1  ,  sina  sin  ?  -  2  (E,  0  +  A*)  sina  sin? 
sinJA  11  lzz 

2  2 

+  E2^  sina  cos  ?  +  E^2  sin  ?  cosa 


-  2E2^  sin?  cos?  cosa+E22  cos^?  cosa] 


'21 


2 


sin3A  ^Ell  s;i-n^  sin  «  +  2  (E12  +  A2)  sina  sin? 


2  2 

+  E21  sin?  cos  a  -  E^2  sin  a  cos? 


-  2E21  sina  cosa  cos?-E:>7  cos  a  cos?] 


22 


"22 


—j-  [E11  sin3a  +  (3E12+2a2)  sin2a  cosa 


sin  A 
?J21 


2  3 

+  3E01  sina  cos  a+E22  cos  a] 


* 

1 

< 

sin? 

-  A* 

U1 

sinA 

L A 

V 

* 

1 

K 

sina 

+  A* 

V 

u2 

sinA 

(F.6) 

(F.  7) 
cos? 

(F .  8) 

cosa 

(F.9) 

(F.10) 

(F.ll) 

(F , 12) 
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APPENDIX  G 

FORMULATION  OF  THE  SYSTEM  OF  EQUATIONS 
FOR  SOLUTION 

The  system  of  equations  (Eq.  (2,25))  can  be 
written  as 


36?  u,  36 

_A  _  -1  _£  =  b 
at  -  s t  i 
q 


(G . la) 


3 

at 


922) 


36 

"at 


* 

2 


96 
_ £ 


(G. lb) 


(G.lc) 


where  ,  b^,  and  b^  are  defined  by  referring  to 

Eq.  (2.25).  Consider  Eq.  (Gl.a)  first.  Using  Eq.  (F.l) 

from  Appendix  F  and  the  definition  of  H.  ,  Eq.  (G.la) 

9p 


becomes 


3  f  1 
3t  sin  A 


* 


sin£ 


cosol 


U1  3 

-  at  (  H0  0n)  "bi 

q  p 


(G .  2) 


Using  Eq.  (E.3)  from  Appendix  E  and  various  shape  factor 
definitions,  some  algebraic  manipulation  gives 


113. 
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sinC 

sinX 


H  + 


cosg 

sinX 


AH  - 


cosg 

sinX 


80 


11 


3t 


cos£ 

sinX 


3Hfi 

_ _£ 

3t 


+  0 


11 


sin£ 

sinX 


cosg 

sinX 


A 


3H 

3t 


cosg 

sinX 


(  H-H0  )  9 


‘  3A  _ 
llj  3t  " 


(G.3) 


Recall  from  Section  2*5  that 


H  =  H(H,  M  ) 


(G.4a) 


such  that 


H0  (H,  M  ) 
P  e 


3H  _  3H  3H  ,  9H_  _^_e 
It  3-  3t  3Mfi  9 t 


3H 


3t 


3H 


6P  3H  + 


3H  8t 


3M 

e 
3 1 


(G.4b) 


(G .  5a) 


(G.5b) 


For  the  case  of  steady  edge  conditions,  3Mg/3t  =  0.  Thus, 
using  the  shape  factor  correlations  given  in  Appendix  C, 
Eqs.  (G.5a,b)  become 


114 


AEDC  TR  83-37 


||  =  (1  +  0.113  M*)  ||  (G, 6a) 

3He 

-Tt2  “  <°-185  Me>  H  (G-6b) 


Substituting  Eqs .  (G.6a,b)  into  Eq.  (G.3)  results  in 


sin£  „  .  cosi 


)  A - —  H 


sinX  sinX 


-  “0 
q  P 


30 


11 


r 


cosg 


sTnX  ASu11'0-072  He>  +|M  9ll,lt0'113  Me> 


3A  . 
3t  bl 


2  U1 

0.185  M“  —  0,, 
6  q  11 

3H  + 
3t 

>®u 

p 

(G.7) 


Analogous  results  can  be  obtained  for  Eqs.  (G.lb,c). 

After  these  operations  have  been  performed f  the  final  form 
of  the  equations  can  be  written  as 


- 

“ 

ail 

a12 

al3 

011 

bi 

a21 

a22 

a23 

3 

3t 

H 

= 

b2 

a31 

a32 

a33 

A 

b3 

(G.8) 


where 
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'll  =  IthI  h  +ffnf  A(H_He  >“~'=  He 

p  q  P 


a  =  A0. . (1-0. 072. M2)  + 5±*4  0  (1+0.113  M2) 

12  sinl  11  e  smA  11  e 


2  U1 

-  0.185  M  —  0,  , 
e  —  11 

q 


13 


Iff  (H-He  )0n 

P 


21 


1  2  2 
[  sin  £  +  sin  a 


•  2. 
sin  A 


+  A(H-Hq  -2)  (cosa  sina  -  cos£  sin£) 

-  A2(H-H0  -1) (cos2£  + cos2a)  ] 

P 


l22  =  - 5“  [  A0  ^  (cosa  sina  -  cos£  sin£) 

sin  A 

-  A201X (cos2£  + cos2a)  ]  (1-0.072  M2) 


23 


— [  0^(H-Hq  -2)  (cosa  sina  -  cos£  sin£) 
sin  A  p 


-  2A0. .  (H-H Q  -1)  (cos2£  +  cos2a)  ] 
11  eo 


a31 


H  -  A(H-H.  )  H. 

smA  sinA  O'  —  9 

P  q  P 


32 


0  .  (1+0.113  M2) 
sinA  11  e 


SIH?  A011(1-O.O72  M2) 


2  U2 

-  0.185  M*  ~  0 

q  11 


(G.9) 

(G.10) 

(G .11) 

(G .  12) 

(G. 13) 

(G.14) 

(G .  15) 

(G . 16 ) 
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'33  -  -  inf  <H-H0p>an 


(G.17) 


bl  * 


-  <J< - ; - ^7 

h1h2  smX  p  q 


L3x1  2 


+  (  hl  SinX  P  q2  °12 


-  cotxf-i  s*  + 


sinX 

- 2 

p  q 

en 

) 

3u1 

•*"  4. 

6* 

3u^ 

hiq 

3x^ 

h2q 

3x2 

+ 

\ 

+ K2  cscX| 


fu. 


12\~  62  +  6 12)  "  2  Cf 
q  /  xi 


(g.18) 


b3 


“  q 


_ 1  r  a 

bih2  sin*  p"  qL^xi 


_  _2 

(  sinX  p  q  ©21  ) 


+  ( hi  sinX  P  q2  ®22  > 


1  i‘  S*  3u2 

rri^rfinq 


u?  *  \  /«!  * 

-  K2  ootx[r|  <!2  +  922j  +  Kx  cscX^—  s1  +  e11/ 


U2  * 

+  k21(-|  61  +  e 


21 


(G.19) 
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b2 


'  2<3\2p  qs  hxh2  sinx  lSil  ‘  "2 

atr  [hl  sinX  P  ^  (  e12  +  £22  >  ]) 


g  _ 3 

[  h0  sinA  p  q  (  +  e21  )  1 


U1  ,  **  .  3ul  1  Su2 

"  ^2  U2  61-U15u2}^ 


U1  /  **  jr*  >  1  ^  1 

,-2  1  u.  ox,  .  -2 
h^q  1  1  n^q 


,  *  *  9ui  u5  *  *  3u2 

+  — L_  (u.  6*  -u,  6*  )  +  (  6*  -6„  ) 

h^-2  1  2  2  ux  73^  h^-2  2  u2  3x2 

__2 

-  K  “tX^+^HJ-6*  )J 

q  J- 

r  u2  i 

-  K2  cot  A  |^e22  +-2  (  6  2  “ 


q 

-2 


+  K 


+  K 


T  U1  *  *  1 

K1  CSCX  |_£12  +^2  (  62  ~  <5U2)J 

-2 

cscA  [e21+z|  K’6u >] 

l  q  i  J 

[£12+^<“l  62*-“2  <J\ 

[e21+=l  in2  6I-“1  Su2>] 


V  U1  “2 

+  (D  +  D  )  )  -  —  b, - b, 

U1  U2  f  q  1  5  3 


U- 


(G .  20) 


(Note  that  b2  and  b3  are  listed  above  in  reverse  order.  This 
is  for  computational  convenience  because  b2  is  defined  in 


terms  of  and  b^) . 
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APPENDIX  H 


LINEAR  STABILITY  ANALYSIS  OF  THE  MODEL  PROBLEM 
USING  FOUR-STAGE  RUNGE-KUTTA  WITH  VARIOUS 
SPACE  DIFFERENCE  APPROXIMATIONS 


First-Order  Backward 


Recall  from  Eq.  (3.13),  the  model  equation  can  be 
written  using  first-order  backward  differences  as 


dUi 

dt 


3.  vi  t  n  \  r  i*  i  n  i 

s  7  (ui 1  -  f(t'ui » 


where 


(H.la) 


n 


n  n 


i  ii  %  n 

V  (  ui  )  =  ui  -  u±_ 


(H.lb) 


Substitution  of  Eq.  (H.lb)  into  Eqs .  (3.3b-e)  results  in 


k-,  =  - 


a  „  ,  n  . 

Ax  V  1  ui  * 


(H.2a) 


.  .  ,  a  At  „2  ,  n  , 

k2  "  kl  +  - 2  7  (  Ui  } 

^  1  2 (Ax) ^  1 


(H.2b) 


k,  = 


k,  -  aiiAu!  v3  <  un  > 

4(Ax)J  1 


(H.2c) 


k  =  -k  +  2k  +  g-  (A^j-  V4  (  u“  ) 
4  1  3  4  { Ax )  4  1 


(H . 2d) 


119 


AEDC-TR-83-37 


where 


2 

V 

,  n 
(ui 

)  = 

n 
u . 
i 

-  2«;_1 

,  n 
+  ui-2 

V3 

(Ui 

)  = 

n 
u  , 
l 

- 3u?-i 

+ 

n 

-  u ,  , 
i-3 

V4 

,  n 
(ui 

)  = 

n 

ui 

- 

+  K-2 

A  11 

-4Ui_ 

(H.3a) 
(H .  3b) 
(H.3c) 


By  letting  the  solution  at  time-level  n  be  represented  by 


where 


n  =  <y(iAx)  _  Ki<f) 
i  e 

=  yAx  (y  =  constant) 
K  =  /T 


(H.4) 


and  use  the  definition  e  +  Kl^  =  cost})  ±  <sin<£,  successive 
substitutions  result  in  the  R-K  scheme  being  expressed  as 


n+1 

Ui 

- —  =  R(CFL,<j>)  +  Kl(CFL,<J>)  (H.5a) 

Ui 

where 

R  ( CFL ,  <p )  =  1  -  (CFL)  (l-costf)  +  \  (CFL)  2  ( l-2cos<j>  +  cos  2  d>) 

-  ^-(CFL)  2  (l-3cos<j>  +  3cos20  -  cos3(fi) 

1  4 

+  2^- (CFL)  (l-4cos<f>  +  6cos2<t>  -  4cos3tj)  +  cos4iJ>)  (H.5b) 
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I(CFL,4>)%  =  -  (CFL)  sin<f>  +  |(CFL)  2  (2sin<j>  -  sin2(f>) 
-  ^ (CFL)  3  (3sin<f>  -  3sin24>  +  sin3cf>) 


1  4 

+  24" (CFL)  (4sin(|)  -  6sin2<f  +  4sin3<f)  -  sin4(})) 


(H.5c) 


CFL  =  a 


(H.5d) 


Thus,  the  amplification  factor  |g[  is 


G  = 


|un+1 

1  l 


[  R2  (CFL  ,  (J>)  +  I2  (CFL ,  (J) )  ]1/2  (H.6) 


which  must  be  less  than  one  for  the  solution  to  remain 
bounded  as  n  -+  00 . 

Similar  expressions  can  be  obtained  for  R(CFL,i|i) 
and  I  (CFL ,  <j>)  for  second-order  central  and  second-order 
backward  spatial  differences,  as  follows. 


Second-Order  Central 


R(CFL,<f>)  -  1  -  ^(CFL)2  sin2<f>  +  ^r(CFL)  ^  sin^tji  (H.7a) 


I tCFL r  (f>)  =  (CFL)  s in<J>  -  |-(CFL)  3  sin 3<J> 


(H.7b) 
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Second-Order  Backward 

R(CFL,4>)  =  1  -  (CFL)  (  |  -  2cos(|)  +|-  cos24>) 

+  (CFL)  2  (  ~  -  3cos<j>  +  A-  cos20  -  cos3<j>  +  y  cos44>) 

-  (CFL)  3  (  A  cosiji  +||  cos 2 tj>  -A  cos3t|> 

19  11 

+  ^  cos44i  -  y  cos5<}>  +  cos64>) 

+  (CFL)  ^  cos#  +  A  cos  2#  -  A  cos  3# 

443  25  9  1 

+  Y92  cos44>  ~Y4  cos5<f>  +  32”  cos6^  ~2T  cos 74» 

+  yp  00s 8 cj)  ) 

I  (CFL,  #)  =  -  (CFL)  (2sin#  -  |  sin2#) 

2  11  1 
+  (CFL)  ( 3  sin#  -  -y  sin2#  +  sin3<|>  -  y  sin4#) 

-  (CFL)  3  (|-  sint})  -  A  sin2tJ)+A  sin3(J>  -  sin4# 

+  i  sin5(}i  -  A  sin6tj)) 

49  81  25  44T 

+  (CFL)  (  y  sintfi  -  jy  sin2i}>  +  sin3#  -  yg-y  sin4<(! 

25  9  1  1 

+  24-  sin5<})  -  jY  Sin6(j>  +  p  sin7<()  -  yp  sin8#  ) 


(H.  8a) 


H.  8b) 
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NOMENCLATURE 


Symbol 


a ,  b 


A 


A 

A 


b 

bl'b2'b3 

c 


Constants  used  in  model  equation, 

Eqs .  (3.1)  and  (3.4) 

Elements  of  A — see  Appendix  G 

Parameter  used  in  Johnston’s  cross-flow 
velocity  profile,  Eq.  (2.31b) 

Coefficient  matrix  of  the  system  of 
equations  to  be  solved,  Eq.  (2.41) — see 
Appendix  G 

Vector  defined  by  Eq.  (2.41)--see 
Appendix  G 

Components  of  b 

Airfoil  chord 


c  Vector  used  in  Eq.  (3.9) 

cf  Skin  friction  coefficient  using  streamwise 

component  of  wall  shear  stress 


Cfx 


i 


CFL 


Skin  friction  coefficient  using  shear 
stress  component  in  x . -direction, 

Eq.  (2.12)  1 

a  At/Ax,  Eq.  (3.17) 


c 

P 

ACFD 


Static  pressure  coefficient 
Parameter  used  in  Appendix  C 


Dx'Dxl'Dx2 


Spatial  difference  operators,  Eqs.  (3.2) 
and  (3. 5) --see  Appendix  H 


DUl rDu?  Dissipation  integrals  expressed  in  the 

nonorthogonal  coordinate  system — 
see  Eq.  (D.l) 

Ds,Dn  Dissipation  integrals,  Eqs.  (H.8c)  and 

(H.8d),  respectively 
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Symbol 
f  ( t,u“? } 


g(u^) 


hl'  h2 


H 

H 


He  * 


H 


0, 


'p 

(Ht) 


rms 


3 


^1' ^2 ' ^3 ' ^4 


K^r  2  #  K2I 


Defined  by  Eq.  (3.2) 

Compressibility  factor,  Eq .  (2.34d) 

Variable  used  in  computing  metric 
coef f icients--see  Appendix  A 

Defined  by  Eq.  (3.5) 

Amplification  factor,  Eq .  (3.16) — see 
Appendix  H 

Surface  metric  coefficients,  Eq.  (A.l) — 
see  Appendix  A 

Streamwise  shape  factor,  Eq.  (2.27b) 

"Incompressible"  shape  factor,  Eq.  (2.27a) 

Streamwise  shape  factor,  Eq.  (2.27d) 

Streamwise  shape  factor,  Eq.  (2.27c) 

Convergence  parameter,  Eq.  (3.19) 

Parameters  corresponding  to  x.,x  - 
directions,  respectively  1  ^ 

Imaginary  part  of  |g|,  Eq.  (3. 16) --see 
Appendix  H 

Velocity  gradient 

Defined  in  Eq.  (3.3)  or  (3.6) 

Surface  curvatures,  Eqs .  (A. 3)  and  A. 4) — 
see  Appendix  A 


L 

L1'L2,L3'L4 


Defined  by  Eq.  (2.21) 

Defined  by  Eq.  (2.22);  also  denotes  body 
length 

Used  in  writing  Eq.  (2.6) — see  Eq.  (2.13) 


m,  n 


Used  as  parameters  in  Appendix  C 
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N,N 


i  j 


NI,NJ 

P 

q 

r 

rw 

R 

Re, 

length 

t 


.Coefficient  matrices  for  3U/3x^  and 
3Uf  -/ vectors,  respectively, 

Eqs.  (3.7)  and  (3.9) 

Coefficient  matrices  for  3U/3x2  and 
3Uj./ax2  vectors,  respectively, 

Eqs.  (3.7)  and  (3.9) 

Maximum  values  of  i  and  j,  respectively 

Static  pressure 

Magnitude  of  total  velocity 

Position  vector — see  Figure  1 

Local  body  radius 

Real  part  of  |g|,  Eq.  (3. 16) --see 
Appendix  H 

Reynolds  number  based  on  some  character¬ 
istic  length;  e.g.,  Re©^ 

Time 


t 


1' 


Defined  by  Eqs.  (D.8a,b) 


At 


T 


L 


u 


i 


uf 


,-ue< 


u  ,u 
s  n 


Time-step  increment 

Defined  in  Eq.  (2.24)  ( i  =  1  or  2) 

Defined  in  Eq.  (2.24) 

Velocity  in  x^-direction;  also  denotes 

general  dependent  variable  in  model 
equation--see  Chapter  III 

Any  component  of  the  vector  of  dependent 
variables  at  the  (i,j)  mesh  point 

Chordwise  and  spanwise  velocity  components, 
respectively,  used  as  inputs  for  Cumpsty 
and  Head  Test  Case--see  Section  5.3 

Velocity  components  in  streamline  coordi¬ 
nate  system 
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Symbol 


—  T 

U  Vector  of  dependent  variables  =  (0  1fH,A) 

--see  Eq.  (2.41)  1 

jj .  .  Vector  of  dependent  variables  at  (i,j) 

mesh  point — see  Eq  .  (3.9) 

xc  Chordwise  coordinate — see  Figures  8  and  9 

x.  General  curvilinear  coordinate 

1  (i  =  1,2,3) — see  Figure  1 

x.  Cartesian  coordinate  (i  =  1,2,3) — see 

1  Figure  1 


Ax 


Either  Ax^  or  Ax2 


Greek  Symbol 

a  Angle  measured  positive  from  xi~axis  to 

local  resultant  edge  velocity  vector 


a 


o 


0 


$ 


w 


Y 


6 


* 

u 


i 


Sweep  angle 

tan  ^{u  /u  ) 
n/  s 

Angle  between  resultant  wall  shear  stress 
vector  and  resultant  edge  velocity  vector 

Ratio  of  specific  heats,  taken  as  1.4; 
also  used  as  an  arbitrary  constant  in 
Appendix  H 

Displacement  thickness  defined  using 
velocities  in  nonorthogonal  coordinate 
system,  Eq.  (2.10a) 

Displacement  thickness  defined  using 
velocities  in  nonorthogonal  coordinate 
system,  Eq.  (2.10d) 

Finite-difference  operator — see  Appendix  H 

Displacement  thickness  defined  using 
velocities  in  streamline  coordinate 
system,  Eqs.  (B.7),  (B.8),  (B.17),  (B.18) 
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Greek  Symbol 


e 


i  j 


E 


i  j 


e.  . 


e 


u 


K 


Energy  thickness  defined  using  velocities 
in  nonorthogonal  coordinate  system, 

Eq.  (2.10c) 

Energy  thickness  defined  using  velocities 
in  streamline  coordinate  system, 

Eqs  .  { B.  13  )-►  ( B .  16  ) 

Momentum  thickness  defined  using  velocities 
in  nonorthogonal  coordinate  systems, 

Eq.  (2.10b) 

Density  thickness,  Eq.  (2.10e) 

Momentum  thickness  defined  using  velocities 
in  streamline  coordinate  system, 

Eqs.  (  B .  9  )-*  ( B .  12  )  ,  (2.27a) 


AT 


A  Angle  between  x,  and  x~  axes,  measured 

positive  from  trie  to  x^  coordinate  axis 

p  Molecular  viscosity 

Pt£  Eddy  viscosity,  Eqs.  (D.12a,b) 

£  X  -  a 

P  Density 

Pf (M) ,  ■  ( N)  Spectral  radii  of  M  and  N  matrices, 

J  J  respectively,  at  (i,j)  mesh  point 

Total  shear  stress  (molecular  plus 
turbulent)  in  x^-direction 

Wave  number 

Smoothing  factor,  Eq.  (3.20) 
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Subscripts 

edge 

i  or  i,j 


s ,  n 


w,  o 


oo 


Denotes  boundary- layer  edge  value 

Denotes  quantity  evaluated  at  the  (i)  or 
(i,j)  mesh  point;  also  denotes  parameters 
corresponding  to  x^-  or  x^-direction 

Denotes  quantity  resolved  in  the  stream¬ 
line  coordinate  system 

Denotes  quantity  evaluated  at  the  body 
surface 

Denotes  quantity  in  the  x^-direction 
Denotes  free-stream  condition 


Superscripts 

(_ ) 


n 


Denotes  boundary-layer  edge  value;  also 
denotes  "incompressible"  quantity 

Denotes  quantity  using  velocity  normal  to 
streamwise  direction 


s 


Denotes  quantity  using  streamwise  velocity 
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